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Abstract:  Model  independent  enhancements  and  adaptations  to  the 
Gauss-Marquardt-Levenberg  (GML)  method  of  computer-based  parame¬ 
ter  estimation  are  described  and  demonstrated  as  potential  improvements 
to  existing  HEC-HMS  automatic  calibration  capabilities.  In  contrast  to  ex¬ 
isting  HEC-HMS  automated  parameter  estimation  capabilities,  these 
methods  support  global  optimization,  the  ability  to  simultaneously  cali¬ 
brate  multiple  subwatershed  systems  represented  within  an  HEC-HMS 
model,  and  they  also  provide  information  about  individual  parameter  sen¬ 
sitivities  and  parameter  correlation  during  and  at  the  end  of  the  calibra¬ 
tion  process.  Moreover,  their  model  independent  nature  allows  one  to  in¬ 
clude  into  the  calibration  process  (l)  state  information  other  than  simply 
stream  discharge  data,  (2)  multiple  periods,  rather  than  a  single  time  win¬ 
dow,  of  the  calibration  dataset(s),  and  (3)  the  ability  to  weight  data  in  or¬ 
der  to  accommodate  a  prediction  specific  calibration  effort  or  to  accom¬ 
modate  suspect  and/or  missing  observations.  The  methods  are 
demonstrated  by  calibrating  HEC-HMS  models  to  subwatershed  systems 
in  the  Goodwin  Creek  Experimental  Watershed. 
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1  Introduction 

Conceptual  model  structures  for  the  continuous  simulation  of  watershed 
hydrology  are  predefined,  prior  to  modeling,  by  the  hydrologist’s  under¬ 
standing  of  the  watershed  system.  With  conceptual  model  structures,  it  is 
not  possible  to  independently  measure  at  least  some  of  the  model  parame¬ 
ters;  hence,  they  must  be  estimated  through  a  formal  model  calibration 
exercise.  Hence,  the  efficacy  of  a  conceptual  model  structure  to  inform  wa¬ 
tershed  management  is  heavily  reliant  upon  observed  system  response 
data  and  the  information  that  one  can  reliably  glean  from  it  during  the 
calibration  process. 

Computer-based  calibration  of  conceptual  watershed  model  structures  for 
continuous  hydrologic  simulation  generally  involves  minimization  of  an 
objective  function  -  a  measure  of  model-to-measurement  misfit.  In  simple 
cases  this  is  comprised  of  differences  between  measured  and  modeled 
flows  at,  for  example,  daily,  hourly,  or  even  smaller  intervals.  In  many 
cases,  observed  and  modeled  flows  are  transformed  (for  example  through 
a  Box-Cox  transformation)  before  fitting,  and/or  residuals  are  fitted  to  an 
ARM  A  model  prior  to  formulation  of  an  objective  function,  in  order  to 
reduce  heteroscedascity  and  temporal  correlation  (Box  and  Tiao  1973,  Box 
and  Jenkins  1976,  Kuczera  1983,  Bates  and  Campbell  2001).  In  more 
complex  cases  a  multi-criteria  objective  function  is  constructed  in  which 
different  measurement  types,  or  the  same  measurement  type  processed  in 
different  ways,  comprises  separate  components  of  a  composite  global 
objective  function  (Madsen  2000,  Boyle  et  al.  2000,  Doherty  and 
Johnston  2003). 

A  unique  solution  to  the  inverse  problem  of  model  calibration  can  only  be 
guaranteed  if  the  number  of  parameters  requiring  estimation  is 
commensurate  with  the  information  content  of  the  calibration  dataset. 
Often  this  is  ensured  by  adherence  to  the  so-called  “principle  of 
parsimony”  in  design  of  the  inverse  problem,  in  which  parameters  requir¬ 
ing  adjustment  through  the  calibration  process  are  reduced  to  a  number 
for  which  a  unique  estimate  can  be  obtained  for  each.  If  calibration  is 
computer-assisted,  then  prior  to  initiating  execution  of  the  parameter  es¬ 
timation  software  package,  the  modeler  selects  those  parameters  that 
he/she  wishes  to  estimate  (normally  on  the  basis  of  anticipated  higher 
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sensitivity  of  model  outputs  to  these  parameters)  and  holds  other  parame¬ 
ters  fixed  at  “sensible  values.” 

Often  this  method  works  well.  However,  problems  associated  with  this  ap¬ 
proach  include  the  following: 

1.  It  is  not  always  possible  to  know  ahead  of  the  parameter  estimation 
process  how  many  parameters  can  be  estimated.  If  too  few  are  selected 
for  estimation  it  may  not  be  possible  to  obtain  a  good  fit  between 
model  outputs  and  field  measurements.  If  too  many  are  selected,  the 
parameter  estimation  process  may  suffer  from  numerical  instability 
and/or  result  in  the  estimation  of  a  set  of  parameter  values  that  lack 
credibility. 

2.  Individual  parameter  sensitivities  are  not  the  sole  determiner  of  what 
is  estimable  and  what  is  not.  Situations  are  often  encountered  where 
model  outputs  have  a  low  sensitivity  to  certain  parameters  collectively, 
but  can  be  very  sensitive  to  the  same  parameters  individually.  This  is 
the  phenomenon  of  “parameter  correlation.” 

3.  Traditional  approaches  to  calibration  are  not  well  suited  to  the  solution 
of  complex  inverse  problems,  such  as  those  involving  simultaneous 
calibration  of  multiple  models,  where  the  estimation  of  useful  values 
for  otherwise  nonuniquely  estimable  parameters  may  be  assisted 
through  the  provision  of  trans-model  parameter  relationships,  from 
which  a  departure  will  only  be  tolerated  if  supported  by  the  calibration 
dataset. 

These  problems  can  be  overcome  through  the  use  of  parameter  estimation 
algorithms  that  allow  mathematical  regularization  to  be  implemented  as 
part  of  the  parameter  estimation  process  itself.  The  result  is  a  stable  solu¬ 
tion  to  the  inverse  problem  (regardless  of  how  ill -posed  it  is),  and  avoid¬ 
ance  of  the  deleterious  effects  of  numerical  instability  on  both  the  parame¬ 
ter  estimation  process  itself,  and  on  the  outcomes  of  that  process,  namely 
the  set  of  estimated  parameter  values.  A  well-designed  regularization  algo¬ 
rithm,  like  its  manual  counterpart,  achieves  numerical  stability  by  re¬ 
formulating  the  inverse  problem  in  a  way  that  recognizes  the  level  of  par¬ 
simony  that  is  necessary  to  attain  a  stable  solution  to  that  problem.  How¬ 
ever,  this  “parsimonizing”  is  undertaken  in  the  context  of  a  specific  cali¬ 
bration  dataset,  allowing  numerical  stability  to  be  achieved  without 
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compromising  model-to-measurement  fit  any  more  than  is  deemed  neces¬ 
sary  by  the  modeler. 

While  measures  can  thus  be  taken  to  ensure  mathematical  tractability  of 
an  inverse  problem  posed  on  the  basis  of  a  properly-processed  calibration 
dataset,  it  is  rarely  possible  to  avoid  the  fact  that  when  calibrating 
conceptual  watershed  model  structures  the  objective  function  will  often 
contain  local  minima  in  addition  to  its  global  minimum  (Wagener,  Whea- 
ter,  and  Gupta  2004,  and  references  cited  therein).  This  presents 
challenges  to  the  design  of  automatic  calibration  software,  for  a  modeler 
who  uses  such  software  has  the  right  to  expect  that  estimated  parameter 
sets  result  in  the  best  possible  fit  between  model  outputs  and  field 
measurements  (with  due  account  taken  of  parameter  believability). 

The  watershed  model  calibration  concepts  noted  above  will  be  explored 
while  pursuing  the  objectives  for  this  article  which  include  describing  and 
demonstrating  the  use  of  parameter  estimation  methodologies  that  could 
potentially  be  employed  to  improve  upon  existing  HEC-HMS  (Hydrologic 
Engineering  Center’s  Hydrologic  Modeling  System)  automated  parameter  es¬ 
timation  capabilities.  In  particular,  the  intent  for  this  article  is  to  describe 
and  demonstrate  the  use  of  methods  that  (1)  accommodate  local  minima 
(Skahill  and  Doherty  2006)  and  (2)  support  the  inversion  of  complex  (i.e., 
highly  parameterized)  models  (Doherty  and  Skahill  2006),  and  through 
this  process  to  also: 

a.  Demonstrate  that  each  of  these  two  methods  provide  information 
about  individual  parameter  sensitivities  and  parameter  correlation. 

b.  Demonstrate  how  state  information  other  than  stream  discharge  data 
alone  can  be  included  into  the  automatic  calibration  process  of  a  HEC- 
HMS  model. 

c.  Demonstrate  how  objective  function  definitions  associated  with  spe¬ 
cific  periods  of  the  calibration  dataset  can  be  included  into  the  auto¬ 
matic  calibration  process  of  a  HEC-HMS  model. 
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2  HEC-HMS  Automatic  Parameter 
Estimation 

HEC-HMS  is  a  widely  used  graphically-based  integrated  watershed  model¬ 
ing  environment  with  a  user  interface  to  a  collection  of  potential  methods 
for  computing  excess  precipitation,  runoff  transformation,  and  hydrologic 
routing,  among  others,  to  support  event-based  or  continuous  simulation  of 
the  hydrologic  response  of  a  dendritic  watershed  system  subject  to  mete¬ 
orological  forcing.  The  graphical  user  interface  enables  seamless  move¬ 
ment  between  the  database,  model  development,  computational  engine(s), 
and  post-processing  capabilities  of  HEC-HMS. 

Existing  HEC-HMS  automatic  calibration  capabilities  include  the  ability  to 
specify  one  of  six  different  objective  functions  and  one  of  two  different 
search  algorithms  (HEC  2005).  Four  of  the  six  objective  function  defini¬ 
tions  currently  available  within  HEC-HMS  are  summarized  in  Table  1 
(HEC  2000). 


Table  1.  Four  of  the  six  objective  function  definitions  available  within  HEC-HMS. 


Description 

Equation 

1-norm 

i1 = Z  l&(0- a  (;)| 

1=1 

2-norm  squared 

*  =  Xle0(i)-eM 

i= 1 

Percent  error  in  peak 

0  =  100 

Qs  {peak)  -  Q0  {peak) 

Qoipeak) 

Peak-weighted  root 
mean  square  error 

_/=i  l  2  Q0{mean)  )  Jj 

CD  =  objective  function;  n  =  number  of  hydrograph  ordinates;  Q0  =  observed  flows;  Qs  =  simu¬ 
lated  flows;  Qo(peak)  =  observed  peak;  Qs(peak)  =  simulated  peak;  and  Q0(mean)  =  mean  of 
observed  flows. 
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The  remaining  two  objective  function  definitions  currently  available 
within  HEC-HMS  include  the  percent  error  in  volume  objective  function 
and  the  time-weighted  objective  function  (HEC  2005).  Within  HEC-HMS, 
the  chosen  objective  function  is  evaluated  for  a  single  specified  period  of 
time. 

The  limited  practical  guidance  provided  in  HEC  (2005)  for  the  selection  of 
a  specific  objective  function  definition  is  supported  by  unrelated  research. 
Conceptual  watershed  model  structural  uncertainty,  due  to  simplifications 
and/or  inadequacies  in  hydrologic  process  descriptions,  often  results  in  a 
model’s  inability  to  fit  all  response  modes  of  the  hydrograph  with  a  unique 
parameter  set  (Wagener,  Wheater,  and  Gupta  2004,  and  references  cited 
therein).  Moreover,  the  predictive  error  variance  analysis  of  Moore  and 
Doherty  (2005)  indicated  that  model  calibration  and  model  prediction 
should  not  be  treated  as  independent  processes. 

The  two  search  algorithms  within  HEC-HMS  include  the  univariate- 
gradient  search  algorithm  and  the  derivative-free  minimization  algorithm 
of  Nelder  and  Mead  (1965). 

The  univariate-gradient  local  search  algorithm  is  implemented  within 
HEC-HMS  in  the  following  manner 


=  x‘  +Ax‘ 

(1) 

(2) 

where  x,  Axk ,  and/represent  the  adjustable  model  parameter,  the  pa¬ 
rameter  upgrade  for  x,  and  the  objective  function,  respectively  (HEC 
2000).  HEC-HMS  approximates  the  derivatives  in  Equation  2  numerically 
using  finite  difference  methods.  If  the  number  of  adjustable  model  pa¬ 
rameters  is  greater  than  one,  this  procedure  is  applied  successively  to  each 
adjustable  parameter  while  holding  all  others  constant. 

The  Nelder  and  Mead  (1965)  local  search  algorithm  is  a  widely  used  de¬ 
rivative-free  minimization  algorithm  that  works  in  multiple  dimensions.  It 
is  sometimes  referred  to  as  an  “amoeba”  method  because  it  works  by  set¬ 
ting  up  rules  that  allow  a  cloud  of  points  in  parameter  space  to  “crawl”  to 
an  objective  function  minimum  in  a  vaguely  amoeboid  fashion.  Rather 
than  starting  with  a  single  initial  guess  for  the  optimized  model,  as  with 
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the  univariate-gradient  search  algorithm,  the  Nelder  and  Mead  method 
begins  by  selecting  m  +  l  parameters  which  form  a  simplex  -  the  simplest 
possible  shape  in  an  m-dimensional  parameter  space.  With  the  Nelder  and 
Mead  simplex  algorithm,  the  m  +  l  vertices  of  a  simplex  of  approximation 
to  an  optimal  point  in  m-dimensional  parameter  space  are  sampled,  or¬ 
dered  by  objective  function  value,  and  an  attempt  made  to  replace  the 
worst  vertex  by  reflection  through  the  convex  hull  of  the  remaining  verti¬ 
ces  using  limited  sampling  along  the  search  direction  so  defined. 

HEC-HMS  supports  constrained  optimization  for  both  search  methods. 
HEC  (2000)  provides  further  details  regarding  the  implementation  of  the 
two  search  algorithms  within  HEC-HMS. 

Some  limitations  associated  with  the  existing  HEC-HMS  automatic  pa¬ 
rameter  estimation  capabilities  include  the  following: 

1.  Existing  search  methods  within  HEC-HMS  are  local;  i.e.,  they  suffer 
from  the  drawback  that  they  may  become  trapped  in  local  objective 
function  minima,  and  thus  report  “optimized”  parameter  values  that 
are  not,  in  fact,  optimized  at  all.  This  can  seriously  degrade  their  utility 
in  the  calibration  of  watershed  models  where  local  optima  abound.  Lo¬ 
cal  search  methods  depend  upon  the  initial  model  estimate.  While  one 
may  attempt  to  rely  on  a  judicious  initial  guess  based  on  physical  rea¬ 
soning  wherein  adjustable  model  parameters  are  derived  from  water¬ 
shed  properties  such  as  soils,  this  approach  suffers  from  the  drawback 
that  the  lumped  conceptual  model  parameters  are  estimated  based  on 
point  samples  analyzed  at  the  laboratory  scale,  and  it  requires  further 
corroboration  (Wagener,  Wheater,  and  Gupta  2004,  and  references 
cited  therein). 

2.  Within  HEC-HMS,  the  calibration  of  multiple  adjacent  gauged  sub¬ 
watersheds  currently  requires  one  to  calibrate  each  subwatershed 
model  independently  of  the  others  rather  than,  for  example,  calibrating 
each  model  individually,  with  due  recognition  of  the  desirability  of  in¬ 
ter- sub  watershed  parameter  similarity  (i.e.,  parameter  values  in  adja¬ 
cent  areas  that  are  associated  with  similar  physiographic  features  rele¬ 
vant  to  hydrologic  response  be  at  least  broadly  similar). 

3.  The  current  methods  within  HEC-HMS  are  not  able  to  report  useful 
information  on  parameter  uncertainty,  correlation,  and  sensitivity  as  a 
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by-product  of  their  use  during  and  after  the  parameter  estimation 
process. 

4.  Objective  functions  can  only  be  comprised  of  differences  between 
measured  and  modeled  flows  (i.e.,  inability  to  construct  a  multi-criteria 
objective  function  in  which  different  measurement  types  (e.g.,  river 
stages,  reservoir  storages,  evapotranspiration,  snowwater  equivalent), 
or  the  same  measurement  type  processed  in  different  ways  (e.g.,  flow, 
baseflow,  quickflow,  volume  aggregations),  comprise  separate 
components  of  a  composite  global  objective  function). 

5.  Objective  functions  can  only  consist  of  a  single  time  window  for  com¬ 
paring  differences  between  measured  and  modeled  flows.  In  addition, 
there  is  no  way  to  weight  data,  say  for  example,  to  guide  a  prediction 
specific  calibration  effort  (Moore  and  Doherty  2005)  or  to  accommo¬ 
date  suspect  and/or  missing  observations. 
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3  Potential  Improvements  for  HEC-HMS 
Automated  Parameter  Estimation 

Approaches 

Enhancements  (Skahill  and  Doherty  2006)  and  adaptations  (Doherty  and 
Skahill  2006)  to  the  Gauss-Marquardt-Levenberg  (GML)  method  of  com¬ 
puter-based  parameter  estimation  (Levenburg  1944,  Marquardt  1963), 
and  a  model  independent  protocol  (Skahill  2006)  wherein  the  inversion 
methods  communicate  with  a  model  through  the  model’s  own  input  and 
output  files,  are  presented  as  approaches  to  address  the  above  noted  limi¬ 
tations  associated  with  existing  HEC-HMS  automatic  parameter  estima¬ 
tion  capabilities. 

Gradient-based  methods  such  as  the  GML  method  have  been  criticized  for 
poor  performance  in  the  face  of  local  optima  (Gupta  et  al.  2003).  Use  of 
such  methods  can  lead  to  the  determination  of  a  parameter  set  that  corre¬ 
sponds  to  a  local,  rather  than  global,  objective  function  minimum,  leaving 
the  user  with  no  idea  of  whether  another  location  exists  within  parameter 
space  for  which  the  objective  function  is  lower.  Nevertheless,  the  method 
also  has  advantages,  chief  among  these  being  its  model-run  efficiency,  and 
its  ability  to  report  useful  information  on  parameter  sensitivities  and  co- 
variances  as  a  by-product  of  its  use.  It  is  also  easily  adapted  to  maintain 
this  efficiency  in  the  face  of  potential  numerical  problems  (that  adversely 
affect  all  parameter  estimation  methodologies)  caused  by  parameter  in¬ 
sensitivity  and/or  parameter  correlation. 

The  present  article  presents  two  algorithmic  enhancements  to  the  GML 
method  that  retain  its  strengths,  but  which  overcome  its  weaknesses  in  the 
face  of  local  optima.  Using  the  first  of  these  methods,  an  “intelligent 
search”  for  better  parameter  sets  is  conducted  in  parameter  subspaces  of 
decreasing  dimensionality  when  progress  of  the  parameter  estimation 
process  is  slowed  either  by  numerical  instability  incurred  through  problem 
ill-posedness,  or  when  a  local  objective  function  minimum  is  encountered. 
The  second  methodology  minimizes  the  chance  of  successive  GML  pa¬ 
rameter  estimation  runs  finding  the  same  objective  function  minimum  by 
starting  successive  runs  at  points  that  are  maximally  removed  from  previ¬ 
ous  parameter  trajectories.  As  well  as  enhancing  the  ability  of  a  GML- 
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based  method  to  find  the  global  objective  function  minimum,  the  latter 
technique  can  also  be  used  to  find  the  locations  of  many  non-global  optima 
(should  they  exist)  in  parameter  space.  This  can  provide  a  useful  means  of 
inquiring  into  the  well-posedness  of  a  parameter  estimation  problem,  and 
for  detecting  the  presence  of  bimodal  parameter  and  predictive  probability 
distributions. 

Moreover,  the  present  article  describes  algorithmic  adaptations  to  the 
GML  method  based  on  an  efficient  and  stable  mathematical  regularization 
scheme.  This  scheme  is  a  variant  of  so-called  “Tikhonov  regularization”  in 
which  the  parameter  estimation  process  is  formulated  as  a  constrained 
minimization  problem.  Use  of  the  methodology  eliminates  the  need  for  a 
modeler  to  formulate  a  parsimonious  inverse  problem  in  which  a  handful 
of  parameters  are  designated  for  estimation  prior  to  initiating  the  calibra¬ 
tion  process.  Instead,  the  level  of  parameter  parsimony  required  to  achieve 
a  stable  solution  to  the  inverse  problem  is  determined  by  the  inversion  al¬ 
gorithm  itself.  Where  parameters,  or  combinations  of  parameters,  cannot 
be  uniquely  estimated,  they  are  provided  with  values,  or  assigned  relation¬ 
ships  with  other  parameters,  that  are  decreed  to  be  realistic  by  the  mod¬ 
eler.  Conversely,  where  the  information  content  of  a  calibration  dataset  is 
sufficient  to  allow  estimates  to  be  made  of  the  values  of  many  parameters, 
the  making  of  such  estimates  is  not  precluded  by  “preemptive  parsimoniz- 
ing”  ahead  of  the  calibration  process. 

Theory 

Gauss-Marquardt-Levenberg  Parameter  Estimation 

Let  the  action  of  a  model  under  calibration  conditions  be  described  by  the 
model  operator  M  that  maps  m-dimensional  parameter  space  to  the  space 
of  the  n  observations  that  are  available  for  use  in  the  calibration  process. 
Let  the  m-dimensional  vector  p  represent  model  parameters  and  the  n- 
dimensional  vector  h  represent  observations.  In  many  instances  of  water¬ 
shed  hydrologic  model  calibration  these  observations  will  represent 
stream  discharges  which  have  been  “processed”  in  some  way  in  order  to 
achieve  homoscedascity,  and  statistical  independence  of  measurement 
“noise.”  The  former  is  often  achieved  through  a  Box-Cox  transformation 
(Box  and  Cox  1964),  while  the  latter  is  often  attempted  through  fitting  re¬ 
siduals  to  an  ARMA  model,  often  as  part  of  the  parameter  estimation 
process  itself  (Box  and  Jenkins  1976,  Kuczera  1983).  The  observations  h 
can  be  comprised  of  a  single  observation  type,  multiple  observation  types, 
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and/or  a  single  observation  type  processed  in  different  ways  in  order  to 
ensure  that  the  information  content  associated  with  different  aspects  of 
the  calibration  dataset  exercise  sufficient  influence  in  the  estimation  of  a 
final  set  of  model  parameters  (Madsen  2000,  Boyle  et  al.  2000,  Doherty 
and  Johnston  2003). 

Model  calibration  seeks  to  minimize  some  measure  of  model-to- 
measurement  misfit  encapsulated  in  a  “measurement  objective  function,” 
herein  designated  as  Om.  In  the  present  instance  this  is  defined  as: 

Om  =  [5W(p)-h]tQ[3f(p)-h]  (3) 

where  Q  is  a  “weight  matrix”  which,  in  the  context  of  watershed  model 
calibration  where  n  is  large,  is  mostly  comprised  of  diagonal  elements 
only.  Ideally,  each  diagonal  element  of  Q  is  proportional  to  the  inverse  of 
the  squared  potential  error  associated  with  the  corresponding  processed 
measurement. 

Where  p  is  estimable  (i.e.,  where  minimization  of  Om  results  in  a  unique 
parameter  set),  it  is  calculated  as: 


p-po  =  (XtQXfiXtQlh-ho)  (4) 

where  X  is  the  model  Jacobian  matrix,  each  row  of  which  is  comprised  of 
the  derivatives  (i.e.,  sensitivities)  of  a  particular  model  output  (for  which 
there  is  a  corresponding  field  measurement)  with  respect  to  all  elements  of 
p.  These  sensitivities  are  calculated  at  current  parameter  values,  repre¬ 
sented  by  p0,  for  which  corresponding  model  outputs  are  h0.  Where  the 
model  is  nonlinear,  p  calculated  through  Equation  4  is  not  optimal  (i.e.,  it 
does  not  minimize  Om)  unless  p0  is  close  to  optimal.  Hence,  after  Equation 
4  is  used  to  calculate  an  improved  parameter  set,  a  new  set  of  sensitivities 
(i.e.,  X)  is  calculated  on  the  basis  of  the  new  parameter  set,  and  the  process 
is  repeated  until  convergence  to  the  objective  function  minimum  is 
achieved. 

In  practice,  the  X‘QX  matrix  of  Equation  4  is  supplemented  by  addition  of 
a  diagonal  term  -  the  so-called  “Marquardt  lambda.”  Thus,  Equation  4  be¬ 
comes: 


p-po  =  (XtQX  +  M)1XtQ(h-h0) 


(5) 
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Normally  X  is  adjusted  during  each  iteration  of  the  parameter  estimation 
process  such  that  its  current  value  results  in  maximum  parameter  im¬ 
provement  during  that  iteration.  When  X  is  high  it  is  easily  shown  that  the 
direction  of  parameter  improvement  is  the  negative  of  the  gradient  of  Om, 
and  under  these  conditions  Equation  5  becomes  equivalent  to  the  “steep¬ 
est  descent”  method  of  parameter  estimation.  While  this  method  can  re¬ 
sult  in  rapid  parameter  improvement  when  parameters  are  far  from  opti¬ 
mal,  its  performance  is  disappointing  in  the  vicinity  of  the  objective 
function  minimum,  especially  where  that  minimum  occupies  a  long  valley 
in  parameter  space  as  a  result  of  excessive  parameter  correlation  or  insen¬ 
sitivity.  In  these  circumstances  “hemstitching”  is  likely  to  occur,  where 
successive  parameter  improvements  result  in  oscillations  across  the  objec¬ 
tive  function  valley,  which  is  never  actually  penetrated.  Hence,  ideally  X 
should  commence  the  parameter  estimation  process  with  a  moderate 
value,  and  then  be  reduced  as  the  process  progresses.  However,  if  X*QX  is 
ill-conditioned,  reducing  the  value  of  X  will  incur  numerical  instability  as 
X*QX  +  AI  of  Equation  5  is  inverted.  Hence,  the  Marquardt  lambda  has  a 
secondary  role,  this  being  that  of  a  de  facto  regularization  device,  with  its 
value  often  being  raised  in  order  to  prevent  instability  in  the  calculation  of 
the  parameter  upgrade  vector  p-p0.  However,  while  the  use  of  a  high 
Marquardt  lambda  can  prevent  a  relatively  ill -posed  parameter  estimation 
problem  from  foundering,  it  achieves  this  at  a  cost  in  efficiency,  for  pa¬ 
rameter  upgrades  become  smaller  at  higher  values  of  X  as  an  inspection  of 
Equation  5  suggests.  Furthermore,  as  stated  above,  the  ability  of  the  cali¬ 
bration  process  to  penetrate  an  elongate  valley  in  parameter  space  may  be 
severely  compromised. 

The  predisposition  of  a  matrix  to  stable  inversion  is  often  measured  by  its 
“condition  number.”  High  condition  numbers  result  in  amplification  of 
numerical  noise  during  the  inversion  process  (Conte  and  de  Boor  1972), 
while  low  condition  numbers  indicate  that  inversion  should  be  possible 
with  little  numerical  difficulty.  In  general,  condition  numbers  for  XlQX 
greater  than  about  ICH  are  to  be  avoided,  for  at  this  level  the  numerical 
noise  incurred  through  finite  difference-based  derivatives  calculation  for 
filling  of  the  X  matrix  is  amplified  to  the  extent  that  parameter  upgrades 
may  lack  integrity.  While  a  raised  Marquardt  lambda  can  often  rescue  such 
a  damaged  process  from  total  failure  as  described  above,  efficiency  of  the 
parameter  estimation  process  is  likely  to  be  seriously  degraded. 
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Another  problem  that  can  be  encountered  when  parameter  estimation  is 
accomplished  by  iterative  calculation  of  p-p0,  using  Equation  5  is  that  this 
process  can  converge  to  a  parameter  set  p  that  corresponds  to  a  local, 
rather  than  the  global,  minimum  of  the  objective  function.  “Gradient 
methods,”  such  as  the  GML  method  described  above,  that  rely  on  equa¬ 
tions  such  as  Equation  5  have  been  criticized  for  this  reason,  and  so-called 
“global  search”  methods  such  as  SCE-UA  (Duan  et  al.  1992)  are  often  used 
instead.  While  a  well-designed  and  robust  global  search  method  can  in¬ 
deed  be  guaranteed  to  minimize  the  objective  function  in  spite  of  the  exis¬ 
tence  of  local  minima,  such  robustness  comes  at  a  price,  this  being  the 
high  number  of  model  runs  that  is  normally  required  for  completion  of  the 
parameter  estimation  process.  To  make  matters  worse,  the  number  of 
model  runs  increases  dramatically  as  the  number  of  parameters  requiring 
estimation  increases.  Use  of  Equation  5,  on  the  other  hand,  is  very  run- 
efficient.  Fortunately,  its  propensity  to  find  local  minima  can  be  mitigated 
through  the  use  of  schemes  such  as  that  described  by  Skahill  and  Doherty 
(2006)  which  combine  the  efficiency  of  gradient  methods  with  the  benefits 
of  introducing  a  small  degree  of  randomness  to  the  parameter  estimation 
process,  together  with  an  ability  to  “learn  from  past  mistakes.”  In  addition, 
Equation  5  can  be  enhanced  by  the  inclusion  of  a  regularization  term 
(much  more  powerful  than  the  Marquardt  lambda  as  will  be  described 
shortly)  that  greatly  increases  the  propensity  for  robust  and  efficient  be¬ 
havior  when  the  dimension  m  of  p  is  large,  and  the  shape  of  the  objective 
function  surface  in  parameter  space  becomes  a  valley  (or  series  of  valleys) 
rather  than  a  bowl  (or  series  of  bowls). 

Gradient-based  methods  such  as  the  GML  method  have  been  criticized  for 
poor  performance  in  the  face  of  local  optima  (Gupta  et  al.  2003).  Use  of 
such  methods  can  lead  to  the  determination  of  a  parameter  set  that 
corresponds  to  a  local,  rather  than  global,  objective  function  minimum, 
leaving  the  user  with  no  idea  of  whether  another  location  exists  within 
parameter  space  for  which  the  objective  function  is  lower.  However  certain 
features  of  the  GML  method  make  it  difficult  to  reject  outright  as  a  serious 
contender  for  use  in  watershed  model  calibration.  These  features  include 
the  following: 

1.  In  calibration  contexts  where  local  optima  are  rare  or  nonexistent,  the 
GML  method  can  normally  find  the  objective  function  minimum  in  far 
fewer  model  runs  than  any  other  method. 
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2.  Estimates  of  parameter  uncertainty,  correlation  and  (in)sensitivity  are 
readily  available  as  a  by-product  of  its  use. 

3.  In  cases  of  high  parameter  insensitivity  and  correlation,  the  method 
can  be  readily  modified  by  the  inclusion  of  various  regularization 
devices  to  maintain  numerical  stability  and  robustness. 

4.  Various  enhancements  can  be  made  to  the  GML  method  that  allow  it  to 
carry  out  linear  or  nonlinear  post-calibration  predictive  uncertainty 
analysis,  with  run  efficiencies  that  far  exceed  those  of  MCMC  methods 
(Vecchia  and  Cooley  1987). 

It  follows  that  if  a  methodology  can  be  found  that  retains  the  advantages  of 
the  GML  method,  while  eradicating  its  propensity  to  be  trapped  in  local 
optima,  such  a  method  would  deserve  serious  consideration  for  use  in 
watershed  model  calibration. 

The  Trajectory  Repulsion  Scheme 

The  robust  performance  of  the  SCE-UA  method,  as  well  as  that  of  most 
other  global  search  methods,  is  based  on  two  principles.  These  are  as 
follows: 

1.  The  injection  of  a  certain  degree  of  randomness  into  the  parameter 
estimation  process  allows  it  to  go  in  directions  that  may  eventually 
prove  fruitful,  even  if  the  attractiveness  of  a  new  direction  may  be 
shielded  by  the  promise  of  local,  more  immediate,  rewards. 

2.  The  benefits  of  randomness  are  partly  offset  by  the  cost  of  making 
mistakes.  Hence,  by  incorporating  into  a  global  optimization  process 
an  ability  to  learn  from  mistakes,  the  likelihood  of  incurring  large  run¬ 
time  penalties  through  repeatedly  making  the  same  (or  a  similar) 
mistake  is  minimized. 

Based  on  these  principles,  a  modified  form  of  the  GML  method  was 
developed  in  order  to  increase  the  capacity  of  this  method  to  work  well  in 
contexts  where  local  minima  occur.  The  package  takes  the  form  of  a  driver, 
in  which  GML  parameter  estimation  is  still  conducted,  but  in  which 
successive  inversion  runs  are  undertaken  under  intelligent  control.  The 
package  is  presently  named  “PD_MS2”  (Skahill  and  Doherty  2006). 
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PD_MS2  commences  execution  by  running  the  model  that  it  must 
calibrate  N  times,  where  N  is  set  by  the  user.  Experience  has  shown  that 
between  the  square  and  the  cube  of  the  number  of  parameters  requiring 
estimation  is  a  suitable  value  for  N.  PD_MS2  employs  random  parameter 
values  for  these  runs;  these  are  sampled  from  a  uniform  or  log-uniform 
distribution  defined  between  user-supplied  upper  and  lower  parameter 
bounds. 

PD_MS2  next  ranks  the  outcomes  of  the  N  random  runs  in  order  of 
increasing  objective  function  value.  It  then  disregards  all  runs  for  which 
the  objective  function  is  above  the  median.  Next,  it  initiates  an  inversion 
run,  with  initial  values  for  this  run  being  equal  to  the  random  parameter 
sample  for  which  the  objective  function  was  lowest.  PD_MS2  monitors 
this  run,  recording  optimized  parameter  values,  as  well  as  parameter 
values  calculated  during  every  iteration  of  the  nonlinear  GML  method 
which  it  implements.  Normally  between  5  and  15  such  iterations  are 
required  to  reach  an  objective  function  minimum.  Each  such  iteration 
requires  that  at  least  as  many  model  runs  be  undertaken  as  there  are 
parameters  requiring  estimation,  plus  a  few  more. 

After  completion  of  the  first  inversion  run,  another  inversion  run  is 
initiated.  For  this  run  it  is  desired  that  the  chances  of  finding  the  same 
objective  function  minimum  as  that  which  was  encountered  on  the  first 
inversion  run  be  minimized.  Hence,  from  among  the  N/ 2  retained  pre¬ 
calibration  samples  of  parameter  space,  a  starting  point  is  chosen  that  is 
maximally  distant  from  any  point  on  the  parameter  trajectory  taken  by  the 
initial  inversion  run.  Selection  of  such  a  starting  point  is  based  on  the 
rationale  that  the  closer  is  a  point  in  parameter  space  to  the  previous 
parameter  trajectory,  the  more  likely  it  is  to  lie  in  the  “catchment  area”  of 
the  previously-encountered  objective  function  minimum. 

After  the  next  inversion  run  is  complete,  another  parameter  set  is  selected 
from  the  N/2  potential  starting  points.  The  parameter  set  selected  is  that 
which  is  maximally  distant  from  all  previous  points  on  all  previous 
trajectories.  The  process  is  then  repeated. 

A  number  of  criteria  can  be  used  to  terminate  the  PD_MS2  global 
optimization  process.  Where  model  run  efficiency  is  an  issue,  PD_MS2 
can  be  instructed  to  cease  execution  if  the  objective  function  has  not  been 
lowered  over  the  last  M,  inversion  runs.  Alternatively,  PD_MS2  can  be 
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asked  to  undertake  M2  inversion  runs  regardless  of  the  outcomes  of  these 
runs.  If  M2  is  moderate  to  large,  this  enables  PD_MS2  to  find  the  locations 
of  many  local  optima  in  parameter  space  (should  these  exist),  thus 
providing  the  user  with  powerful  insights  into  the  structure  of  the 
objective  function  surface. 

It  is  worth  noting  that,  as  well  as  providing  insights  into  the  “broadscale” 
structure  of  the  objective  function  response  surface,  PD_MS2  provides 
insights  into  the  structure  of  this  surface  in  the  vicinity  of  the  global 
objective  function  minimum  as  well.  As  has  already  been  mentioned,  the 
GML  method  can  provide  parameter  sensitivities  and  can  calculate  a  linear 
approximation  to  the  parameter  covariance  matrix,  as  well  as  statistics 
derived  from  this  matrix  including  correlation  coefficients  and 
eigenvectors/eigenvalues  of  the  covariance  matrix.  Information  of  this 
type  is  forthcoming  only  with  difficulty  from  global  search  methods,  this 
difficulty  increasing  with  the  number  of  parameters  being  estimated  and 
with  the  degree  of  correlation  between  them  (which,  unfortunately,  is  the 
very  situation  in  which  such  information  is  of  most  value). 

Temporary  Parameter  Immobilization 

“Temporary  parameter  immobilization”  can  be  used  as  both  a  regulariza¬ 
tion  device  and  as  a  device  for  conducting  ordered  attempts  to  break  out  of 
local  pits  in  parameter  space.  This  scheme  is  implemented  only  if  the  ob¬ 
jective  function  improvement  attained  during  a  particular  iteration  of  the 
GML  process  is  less  than  a  user-supplied  threshold  (normally  to  percent). 
In  implementing  this  scheme,  the  most  insensitive  parameter  is  selected, 
and  temporarily  removed  from  the  optimization  process.  With  the  dimen¬ 
sionality  of  estimable  parameter  space  thus  reduced  (and  with  the  most 
troublesome  parameter  being  temporarily  removed  from  the  parameter 
estimation  process),  the  parameter  upgrade  vector  (which  now  has  no 
component  in  the  subspace  of  parameter  space  occupied  by  the  temporar¬ 
ily  frozen  parameter)  is  recalculated  using  Equation  5.  A  model  run  is  then 
conducted  on  the  basis  of  the  trial  parameter  set  thus  calculated  in  order 
to  compute  the  objective  function  associated  with  this  parameter  set. 
Unless  the  objective  function  has  fallen  by  a  significant  amount,  the  next 
most  troublesome  parameter  is  temporarily  frozen  (in  addition  to  the 
first),  and  the  parameter  upgrade  calculation  procedure  is  repeated.  After 
a  number  of  parameters  have  been  successively  frozen  in  this  manner 
(with  already  frozen  parameters  maintained  in  their  frozen  state),  the 
process  is  abandoned,  and  then  recommenced  using  a  different  value  of 
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the  Marquardt  lambda.  For  a  parameter  estimation  problem  involving  m 
parameters,  up  to  half  of  these  parameters  may  be  progressively  frozen  for 
up  to  three  Marquardt  lambdas,  this  requiring  3777/2  model  runs  for  that 
iteration  for  the  testing  of  parameter  upgrade  vectors  in  addition  to  the 
(depending  on  whether  forward  differences  or  central  differences  are  em¬ 
ployed)  m  or  2 m  model  runs  required  for  filling  of  the  Jacobian  matrix. 
(Note,  however,  that  the  process  is  immediately  abandoned  if  a  suitable 
objective  function  improvement  is  obtained.)  Thus,  implementation  of  the 
TPI  process  may  lead  to  the  requirement  that  between  twice  and  three 
times  (at  the  very  most)  the  number  of  model  runs  be  carried  out  com¬ 
pared  to  normal  GML  operations.  However,  experience  has  demonstrated 
that  on  most  occasions  in  which  the  TPI  method  is  employed  about  50  per¬ 
cent  extra  model  runs  need  to  be  carried  out,  and  that  this  is  generally  a 
small  price  to  pay  for  the  benefits  that  it  brings  in  terms  of  increased  nu¬ 
merical  stability  in  situations  of  parameter  nonuniqueness,  and  for  a  dra¬ 
matic  reduction  in  the  risk  of  becoming  trapped  in  local  objective  function 
pits. 

The  decreased  probability  of  ensnarement  in  local  optima  that  attends  use 
of  the  TPI  scheme  has  its  roots  in  a  number  of  properties  of  this  scheme. 
One  obvious  reason  for  a  heightened  probability  of  success  in  finding  its 
way  out  of  small  regions  of  attraction  of  limited  extent  in  parameter  space 
is  the  sheer  number  of  parameter  upgrades  that  are  attempted  by  this 
scheme,  together  with  the  fact  that  the  directions  pertaining  to  these  up¬ 
grade  attempts  tend  to  be  maximally  different  with  respect  to  each  other. 
This  maximality  of  difference  is  a  result  of  two  factors.  The  first  is  the  fact 
that  the  upgrade  direction  tends  to  be  dominated  by  insensitive  parame¬ 
ters  where  all  parameters  are  involved  in  the  computation  of  this  direc¬ 
tion;  this  is  a  direct  result  of  the  fact  that,  because  of  their  insensitivity,  the 
GML  parameter  estimation  algorithm  calculates  that  these  parameters  re¬ 
quire  larger  movement  than  other  parameters  to  affect  the  objective  func¬ 
tion.  As  dimensions  of  parameter  space  are  progressively  closed  to  the  pa¬ 
rameter  upgrade  vector  through  the  temporary  immobilization  of 
insensitive  parameters,  and  new  upgrade  directions  are  accordingly  com¬ 
puted  in  spaces  of  lower  dimensions,  these  new  directions  will  tend  to  be 
orthogonal  to  the  original  upgrade  vector  which  was  dominated  by  the 
now-omitted  dimensions.  The  penchant  for  orthogonality  is  further  in¬ 
creased  as  a  result  of  the  fact  that  the  entire  dimensionality  reduction 
process  is  repeated  for  widely  different  Marquardt  lambda  values.  As 
documented  in  works  such  as  Bard  (1974),  computed  upgrade  directions 
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can  vary  between  that  of  steepest  descent  down  the  objective  function  sur¬ 
face  when  the  Marquardt  lambda  is  high,  to  a  direction  that  can  be  almost 
orthogonal  to  this  when  the  Marquardt  lambda  is  low. 

Another  important  factor  behind  the  success  of  the  TPI  scheme  is  that  it 
lowers  the  chances  of  upgraded  parameters  finding  local  optima  in  the 
first  place.  Unless  objective  function  improvement  during  a  particular  it¬ 
eration  is  acceptably  large  without  the  help  of  the  TPI  scheme  (which  often 
occurs  in  the  early  stages  of  the  parameter  estimation  process),  use  of  the 
TPI  scheme  requires  that  model  runs  be  carried  out  specifically  to  test  the 
ability  of  different  upgrade  vectors  (often  with  very  different  directions  as 
discussed  above)  to  lower  the  objective  function.  The  upgrade  vector  that 
results  in  the  largest  objective  function  decline  is  that  which  is  selected  as 
the  basis  for  the  next  linearization  of  the  inverse  problem.  Of  all  the  up¬ 
grade  vectors  tested,  this  is  the  one  least  likely  to  lead  to  a  local  objective 
function  minimum,  for  the  encroachment  of  global  or  local  optimality  (for 
which  derivatives  of  the  objective  function  with  respect  to  all  model  pa¬ 
rameters  is  zero)  is  normally  marked  by  smaller  and  smaller  declines  in 
the  objective  function  per  iteration  as  the  GML  method  ensures  that  a  pa¬ 
rameter  set  is  found  from  which  all  directions  lead  uphill.  In  fact,  the  more 
nonlinear  is  the  problem,  the  less  likely  it  is  that  a  parameter  upgrade  vec¬ 
tor  resulting  in  a  large  objective  function  decline  will  lead  directly  to  the 
bottom  of  an  objective  function  minimum  (due  to  the  fact  that  the  equa¬ 
tions  upon  which  this  upgrade  vector  are  calculated  are  based  on  an  as¬ 
sumption  whose  inapplicability  grows  with  increasing  parameter  move¬ 
ment,  and/or  increasing  changes  in  model  outputs  on  account  of  this 
movement). 

An  additional  factor  that  contributes  to  the  success  of  the  TPI  scheme  in 
both  avoidance  of  local  minima  of  small  lateral  extent,  and  in  extricating 
itself  from  such  minima,  is  use  of  finite  differences  for  parameter  deriva¬ 
tives  calculation.  As  was  mentioned  above,  parameter  increments  of  one 
percent  are  often  employed  for  forward  difference  derivatives  calculation 
and  two  percent  for  central  difference  derivatives  calculation.  These  in¬ 
crements  are  large  enough  to  “see”  outside  of  a  small  pit  in  which  it  may  be 
currently  trapped.  Alternatively,  if  current  parameter  values  lie  just  out¬ 
side  of  a  small  pit,  these  increments  are  large  enough  for  the  effect  of  the 
pit  to  exert  a  smaller  influence  on  calculated  derivatives  than  would  be  the 
case  if  derivatives  were  exact.  Thus,  the  use  of  finite-difference-based  pa¬ 
rameter  derivatives  provides  a  kind  of  filtering  mechanism  through  which 
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finer  details  of  the  objective  function  surface  are  prevented  from  conceal¬ 
ing  the  broader  features  of  that  surface. 

So,  through  a  combination  of  the  fact  that  many  upgrade  vectors  are 
tested,  that  a  parameter  upgrade  selection  procedure  is  adopted  that 
minimizes  the  chances  of  being  trapped  in  a  local  minimum  in  the  first 
place,  and  maximizes  the  chances  of  escaping  from  that  minimum  if  en- 
snarement  does  indeed  occur,  and  because  parameter  upgrades  possess 
some  immunity  to  the  effects  of  pits  because  their  calculation  is  based  on 
finite-difference  derivatives  rather  than  point  derivatives,  use  of  the  TPI 
method  in  calibration  of  surface  water  models  has  consistently  resulted  in 
good  performance  in  estimating  parameters  for  those  models. 

(Note  that  selection  of  a  TPI  activation  threshold  of  to  percent  improve¬ 
ment  in  the  objective  function  is  somewhat  arbitrary.  However,  experience 
has  demonstrated  that  this  normally  results  in  efficient  implementation  of 
the  method.  If  the  threshold  is  set  too  high,  TPI-based  parameter  upgrade 
recomputation  will  be  undertaken  on  most  GML  optimization  iterations, 
irrespective  of  proximity,  or  otherwise,  to  an  objective  function  minimum. 
This  can  result  in  wasted  model  runs  if  rapid  objective  function  improve¬ 
ment  is  taking  place  without  the  need  for  TPI  upgrade  repetitions.  On  the 
other  hand,  if  the  improvement  threshold  is  set  too  low,  then  needless 
“struggling”  of  the  GML  method  in  the  face  of  difficulties  incurred  through 
problem  ill-posedness  or  proximity  to  a  local  minimum,  resulting  in  only 
small  improvements  in  the  objective  function  in  successive  iterations,  can 
be  avoided.) 

Regularized  Inversion 

Conceptually,  singularity  or  near-singularity  of  XlQX  (as  occurs  when 
large  numbers  of  parameters  require  estimation  and/or  when  the  informa¬ 
tion  content  of  the  calibration  dataset  with  respect  to  estimated  parame¬ 
ters  is  poor)  can  be  remedied  through  the  addition  of  extra  “observations” 
to  the  parameter  estimation  process  which  pertain  directly  to  the  parame¬ 
ters  requiring  estimation.  For  example,  it  may  be  “observed”  that  each  pa¬ 
rameter  is  equal  to  a  certain,  user-supplied  value;  presumably  this  value 
will  have  been  chosen  to  be  realistic  in  terms  of  the  system  property  which 
the  parameter  represents.  Alternatively  (or  as  well),  it  may  be  “observed” 
that  certain  pairs  of  parameters  are  equal,  or  have  values  which  observe  a 
certain  ratio  or  difference. 
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Let  these  “regularization  relationships”  be  represented  by  the  operator  Z 
acting  on  the  parameter  set  p,  and  let  the  “observed”  values  of  these  rela¬ 
tionships  be  represented  by  j.  Then  the  regularization  relationships  (also 
referred  to  as  “regularization  constraints”  herein)  can  be  represented  by 
the  equation: 


Z{  p)=j 


(6a) 


the  linearized  form  of  which  is: 


Zp  =  j  (6b) 

where  Z  is  the  Jacobian  of  the  Z  operator.  Note  that,  as  is  discussed  below, 
it  is  not  essential  that  Equations  6a  and  6b  be  exactly  observed,  only  that 
they  be  observed  to  the  maximum  extent  possible  in  calibrating  the  model. 

If  the  regularization  constraints  are  given  sufficient  weight  in  comparison 
with  the  observation  weights  encapsulated  in  Q,  a  well-posed  inverse  prob¬ 
lem  will  have  been  formulated.  Mathematically,  this  problem  is  then  itera¬ 
tively  solved  for  the  parameters  p  using  the  equation: 

p-po  =  (XtQX+  p2ZtSZ  +  M)-1(XtQ[h-h0]+  p^SU-jo])  (7) 

In  Equation  7,  j0  represents  the  right  side  of  Equation  6a  when  current  pa¬ 
rameter  values  p0  are  substituted  for  p  in  this  equation.  S  is  a  “relative 
weight  matrix”  assigned  to  the  regularization  observations  j;  it  has  the 
same  role  for  regularization  observations  as  Q  does  for  field  observations. 
All  of  the  relative  regularization  weights  encapsulated  in  S  are  multiplied 
by  a  “regularization  weight  factor”  P2  in  Equation  7  prior  to  calculation  of 
P-Po. 


Selection  of  an  appropriate  value  for  P2  is  critical.  If  its  value  is  too  high 
the  parameter  estimation  process  will  ignore  the  measurement  dataset  h 
in  favor  of  fitting  the  regularization  observations  j.  If  it  is  too  small,  the 
regularization  observations  will  not  endow  the  parameter  estimation  proc¬ 
ess  with  the  numerical  stability  which  it  needs  in  order  to  obtain  estimates 
for  the  parameters  p. 

Equation  7  can  be  shown  to  constitute  a  constrained  minimization  prob¬ 
lem  in  which  a  “regularization  objective  function”  O,  defined  as: 
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Or=[Z(p)-j]tS[Z(p)-j]  (8) 

is  minimized  subject  to  the  constraint  that  Om  of  Equation  3  rises  no 
higher  than  a  user-specified  value,  referred  to  herein  as  the  “target  meas¬ 
urement  objective  function.”  Thus  the  user  informs  the  regularized  inver¬ 
sion  process  of  the  level  of  model-to-measurement  misfit  required;  this 
process  then  enforces  the  regularization  constraints  defined  through 
Equation  6a  to  the  maximum  extent  that  it  can  by  minimizing  O,  subject  to 
the  constraint  that  Om  rises  no  higher  than  the  target  level.  If  the  target 
measurement  objective  function  cannot  be  achieved,  the  regularized  inver¬ 
sion  process  simply  minimizes  Om;  however,  where  minimization  of  Om 
would  otherwise  be  an  unstable  process  due  to  parameter  nonuniqueness, 
stability  of  this  process  is  maintained  by  seeking  that  set  of  parameters  ly¬ 
ing  within  the  elongate  Om  valley  that  also  minimizes  Or.  In  either  case, 
the  regularization  weight  factor  (32  can  be  viewed  as  a  Lagrange  multiplier 
associated  with  the  constrained  minimization  problem,  and  it  is  recalcu¬ 
lated  during  every  iteration  of  the  regularized  nonlinear  parameter  estima¬ 
tion  process  using  a  bisection  algorithm  based  on  local  linearization  of  the 
constrained  minimization  problem  about  current  parameter  values. 

Note  the  continued  inclusion  of  the  Marquardt  lambda  in  Equation  7.  Its 
value  is  adjusted  as  needed  from  iteration  to  iteration  as  a  practical  meas¬ 
ure  to  enhance  optimization  efficiency  and  to  ensure  stability  of  the  pa¬ 
rameter  estimation  process  should  XlQX+  (32ZtSZ  become  ill-conditioned 
through  use  of  an  inappropriately  low  value  for  (32.  This  can  occur  where 
regularization  constraints  are  poorly  formulated,  or  where  too  good  a  fit  is 
sought  between  model  outputs  and  field  measurements,  requiring  that 
regularization  constraints  be  abandoned  in  pursuit  of  this  fit.  Often  it  oc¬ 
curs  for  a  combination  of  these  reasons,  where  weights  on  some  regulari¬ 
zation  constraints  must  be  lowered  for  attainment  of  a  good  fit  between 
model  outputs  and  field  measurements,  but  where  the  relaxation  of  regu¬ 
larization  constraints  then  leads  to  unestimability  of  those  model  parame¬ 
ters  whose  estimation  is  not  realized  through  attainment  of  this  fit. 

Formulation  of  the  inverse  problem  as  a  constrained  minimization  prob¬ 
lem  through  use  of  Equation  7  allows  many  more  parameters  to  be  esti¬ 
mated  than  would  otherwise  be  possible,  thereby  ensuring  that  maximum 
information  is  extracted  from  the  calibration  dataset.  If  the  relationships 
of  Equation  6  are  realistic,  the  fact  that  estimated  parameters  are  such  as 
to  ensure  minimal  deviation  from  these  relationships  heightens  the  prob- 
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ability  that  estimated  parameters  will  themselves  be  realistic.  However,  a 
practical  problem  that  is  often  encountered  when  using  the  Tikhonov 
method  is  that  the  regularization  weight  matrix  S  must  be  supplied  ahead 
of  the  regularized  inversion  process;  furthermore,  it  is  not  adjusted 
through  this  process  except  for  global  multiplication  by  (32.  Ideally,  indi¬ 
vidual  regularization  constraints  described  by  the  rows  of  Equation  6 
should  be  more  strongly  enforced  where  the  information  content  of  the 
calibration  dataset  is  insufficient  to  require  their  contravention  for  the 
sake  of  obtaining  an  appropriate  level  of  model-to-measurement  fit.  How¬ 
ever  because  it  is  almost  impossible  to  know  ahead  of  the  calibration  proc¬ 
ess  the  extent  to  which  this  should  occur  for  each  of  the  different  relation¬ 
ships  encapsulated  in  Z,  it  is  often  very  difficult  to  supply  an  S  matrix  that 
is  an  appropriate  complement  to  the  current  calibration  dataset. 

Adaptive  Regularization 

An  “adaptive  regularization”  methodology  is  now  presented  which  over¬ 
comes  this  problem  in  many  modeling  contexts.  The  set  of  regularization 
constraints  described  by  Equation  6  is  subdivided  into  groups;  if  desired, 
each  constraint  can  be  assigned  to  its  own  group.  The  set  of  model  pa¬ 
rameters  p  is  then  supplemented  by  an  additional  parameter  set  pr,  with 
one  new  parameter  being  defined  for  each  new  regularization  group.  Each 
such  parameter  is,  in  fact,  the  inverse  of  a  group-specific  regularization 
weight  multiplier;  this  group-specific  weight  multiplier  is  applied  in  addi¬ 
tion  to  the  global  weight  multiplier  (32  depicted  in  Equation  7,  the  latter  be¬ 
ing  adjusted  as  part  of  the  constrained  minimization  process  as  described 
above.  Regularization  constraints  are  then  provided  for  the  elements  of  pr  - 
so  that  these  too  can  be  estimated  as  part  of  the  regularized  inversion 
process.  Each  such  constraint  comprises  the  “observation”  that  the  respec¬ 
tive  element  of  pr  is  zero. 

The  reformulated  regularized  inversion  problem  remains  a  constrained 
minimization  process,  and  thus  still  seeks  to  find  a  parameter  set  that  ei¬ 
ther  minimizes  the  measurement  objective  function  Om,  or  reduces  it  to  a 
user-specified  target  level,  while  ensuring  that  the  regularization  objective 
function  Or  is  conditionally  minimized.  Because  conditional  minimization 
of  the  regularization  objective  function  now  requires  maximization  of 
weights  assigned  to  individual  or  groups  of  regularization  constraints, 
these  weights  are  applied  as  strongly  as  possible,  thereby  maximizing  the 
extent  to  which  the  corresponding  regularization  relationships  encapsu¬ 
lated  in  Equation  6  are  adhered.  However,  with  the  calculation  of  the  over- 
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all  regularization  weight  factor  (32  by  the  constrained  minimization  process 
being  such  as  to  allow  minimization  of  the  target  measurement  objective 
function,  or  achievement  of  a  user-specified  target  for  this  function,  these 
regularization  constraints  are  not  so  strongly  enforced  that  model-to- 
measurement  fit  is  compromised.  Thus,  the  regularized  inversion  process 
itself  ensures  that  the  strength  of  enforcement  of  regularization  con¬ 
straints  on  parameter  values  or  relationships  complements  the  informa¬ 
tion  content  of  the  calibration  dataset  in  relation  to  these  parameters.  As  a 
result,  regularization  constraints  are  automatically  applied  more  strongly 
where  the  attainment  of  a  satisfactory  level  of  model -to-measurement  fit 
does  not  require  otherwise,  thus  overcoming  a  disadvantage  of  the  Tik¬ 
honov  method.  The  outcome  is  a  numerically  stable  regularized  inversion 
process  that  achieves  a  desired  level  of  model-to-measurement  fit  with 
impressive  run  economy,  and  that  yields  sensible  values  for  model  pa¬ 
rameters. 

Like  all  numerical  strategies,  this  adaptive  regularization  methodology  is 
more  suitable  for  use  in  some  contexts  than  in  others.  It  is  certainly  not 
the  only  means  by  which  numerical  stability  of  a  regularized  inversion 
process  can  be  achieved,  for  so-called  “subspace  methods”  (Aster  et  al. 
2005)  are  very  effective  in  this  regard.  However,  use  of  the  present  meth¬ 
odology  can  be  beneficial  in  those  modeling  contexts  where  the  means  by 
which  numerical  stability  is  achieved  is  just  as  important  as  the  achieve¬ 
ment  of  that  stability  itself.  In  general,  where  the  necessity  for  parameters 
to  observe  key  values  or  relationships  to  the  maximum  extent  possible 
without  compromising  fit  between  model  outputs  and  field  measurements 
is  a  critical  part  of  the  calibration  process,  then  the  adaptive  regularization 
methodology  described  herein  will  serve  that  calibration  process  well;  such 
a  case  is  demonstrated  in  the  following  section.  However,  the  need  to  in¬ 
troduce  extra  parameters  into  the  calibration  process  in  order  to  guarantee 
enforcement  of  desired  parameter  relationships  does  place  some  restric¬ 
tions  on  the  method.  Where  such  relationships  fall  into  a  relatively  small 
number  of  distinct  groups,  and/or  where  the  number  of  parameters  re¬ 
quiring  estimation  is  not  such  as  to  introduce  vastly  different  levels  of  “es- 
timability”  between  them  (thus  requiring  the  introduction  of  many  new 
parameters  in  order  to  accommodate  the  differential  strengths  with  which 
regularization  constraints  must  be  applied),  the  above  method  has  proven 
very  successful.  However,  where  large  numbers  of  parameters  require  es¬ 
timation,  and  where  differences  in  estimability  between  them  are  likely  to 
cover  a  broad  range,  recourse  to  subspace  methods  becomes  a  necessity. 
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Unfortunately,  in  this  case,  the  guarantee  of  numerical  stability  that  ac¬ 
companies  use  of  such  methods  is  attained  at  the  cost  of  loss  of  ability  on 
the  part  of  the  modeler  to  insist  on  the  observance  of  specified  parameter 
relationships  in  attaining  that  stability. 
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4  Examples 

Use  of  the  methodologies  discussed  in  the  preceding  section  are  now  dem¬ 
onstrated  by  applying  them  to  the  calibration  of  two  different  HEC-HMS 
models  deployed  to  the  Goodwin  Creek  Experimental  Watershed.  Good¬ 
win  Creek  is  an  8.26-square-mile  experimental  watershed  (upland  erosion, 
instream  sediment  transport,  and  watershed  hydrology)  operated  by  the 
Agricultural  Research  Service  (ARS)  of  the  United  States  Department  of 
Agriculture  (USDA).  The  Goodwin  Creek  watershed  is  located  in  northern 
Mississippi,  approximately  60  miles  south  of  Memphis,  TN.  Goodwin 
Creek  is  divided  into  14  nested  subwatersheds  with  a  flow  measuring 
flume  constructed  at  each  of  the  subwatershed  outlets.  The  drainage  areas 
above  the  gauging  sites  range  from  0.63  to  8.26  square  miles.  Thirty-one 
standard  recording  rain  gauges  are  uniformly  located  within  and  just  out¬ 
side  of  the  watershed.  For  further  details  about  the  Goodwin  Creek  Ex¬ 
perimental  Watershed,  see  Blackmarr  (1995).  Digital  elevation  model- 
derived  subwatershed  boundaries  and  stream  network,  and  rain  and 
stream  gauge  locations  are  shown  in  Figure  1. 

The  first  HEC-HMS  model  was  applied  to  the  39.8-acre  drainage  area 
above  streamflow  gauging  station  number  9  in  Goodwin  Creek.  Precipita¬ 
tion  data  from  two  nearby  gauges  and  evaporation  data  associated  with 
four  surrounding  locations  were  supplied  as  meteorological  forcing  data  to 
support  HEC-HMS  continuous  hydrologic  simulation  using  the  deficit 
constant  loss  method;  for  details  about  the  deficit  constant  loss  method 
within  HEC-HMS,  see  HEC  (2005). 


ERDC/CHL  TR-06-13 


25 


Figure  1.  Goodwin  Creek  Experimental  Watershed,  delineated  subwatersheds,  derived 
stream  network,  rain  and  numbered  streamflow  gauge  locations. 


Estimation  of  six  HEC-HMS  parameters  was  undertaken  by  matching  ob¬ 
served  and  simulated  flow  data  over  17  non-contiguous  time  intervals 
spanning  the  period  1  Jan  1989  to  29  Feb  1992.  The  17  periods  were  identi¬ 
fied  based  on  a  manual  inspection  of  the  observed  flow  data.  The  calibra¬ 
tion  period  1  Jan  1989  to  29  Feb  1992  was  selected  based  on  the  determi¬ 
nation  that  there  was  no  land  use/land  cover  alteration  within  the 
subwatershed  for  the  period  1987  through  1992.  Two  calibration  experi¬ 
ments  were  performed  with  the  first  HEC-HMS  model.  The  first  calibra¬ 
tion  experiment  involved  comparing  only  simulated  and  observed  flows; 
the  second  calibration  experiment  involved  comparing  simulated  and  ob¬ 
served  flows  and  base  flows.  Over  the  17  non-contiguous  time  intervals 
spanning  the  calibration  period  1  Jan  1989  to  29  Feb  1992,  the  first  ex¬ 
periment  involved  matching  simulated  flows  with  observed  flows  above  a 
predetermined  threshold  (which  was  uniquely  determined  for  each  of  the 
17  time  intervals  based  on  a  manual  inspection  of  the  data);  whereas,  the 
second  experiment  involved  matching  not  only  simulated  flows  with  ob¬ 
served  flows  (above  the  predetermined  threshold),  but  also  matching 
simulated  base  flows  with  observed  base  flows  for  the  same  times  simu¬ 
lated  and  observed  flows  were  matched.  Hence,  the  first  experiment  re¬ 
sulted  in  a  total  of  2,053  flow  observations  for  use  in  the  calibration  proc- 
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ess;  whereas,  the  second  experiment  resulting  in  a  total  of  4,106  observa¬ 
tions  for  use  in  the  calibration  process. 

In  each  case,  the  objective  function  was  defined  as  the  sum  of  weighted 
squared  differences  between  modeled  and  observed  log-transformed  flow 
data,  with  all  weights  assigned  a  value  of  1.0.  Thus  h  of  Equation  3  was 
comprised  of  the  logs  of  flow  data,  while  the  model  represented  by  X  in 
these  equations  calculated  the  model-generated  counterparts  to  these 
logged  flow  data.  Q  was  the  identity  matrix.  For  both  experiments  with  the 
first  HEC-HMS  model,  the  simulation  time  interval  was  one  hour,  which 
equaled  the  temporal  resolution  of  the  input  precipitation  data. 

Table  2  lists  the  17  non-contiguous  time  intervals  and  the  specified  thresh¬ 
olds.  Figure  2  is  a  plot  of  the  observed  flow  data  that  were  compared  with 
simulated  flows  to  calibrate,  for  both  experiments,  the  first  HEC-HMS 
model.  The  “observed”  base  flow  data  for  comparing  with  HEC-HMS 
simulated  base  flow,  for  the  second  experiment  with  the  first  model,  was 
determined  using  a  quickflow  digital  filter  (Nathan  and  McMahon  1990), 
with  a  value  of  0.999  specified  for  the  scaling  parameter.  For  purposes  of 
illustration,  base  flow  “observations,”  together  with  the  “observed”  quick- 
flow  and  the  total  observed  flow,  for  a  specific  period  during  the  first  of  the 
17  time  intervals,  are  plotted  in  Figure  3.  Examining  Figure  3,  clearly  the 
perceptual  model  (Beven  2001)  was  formulated  to  be  consistent  with  past 
research  (Downer  and  Ogden  2003)  in  that  the  watershed  system  is  domi¬ 
nated  by  quickflow  response. 


Table  2.  Non-contiguous  time  intervals  and  thresholds,  both  which  determined  model  to 
measurement  misfit  for  both  experiments  with  the  first  model. 


Number 

Time  Interval 

Threshold  (cfs) 

1 

01/02/1989  00:00:00  -  02/13/1989  23:00:00 

0.1 

2 

02/18/1989  06:00:00  -  02/21/1989  12:00:00 

0.1 

3 

02/27/1989  02:00:00  -  05/17/1989  23:00:00 

0.1 

4 

05/22/1989  00:00:00  -  08/08/1989  23:00:00 

0.1 

5 

01/03/1990  18:00:00  -  01/08/1990  08:00:00 

0.1 

6 

01/28/1990  20:00:00  -  01/29/1990  23:00:00 

0.1 

7 

02/09/1990  13:00:00  -  03/13/1990  23:00:00 

0.1 

8 

04/06/1990  00:00:00  -  04/06/1990  23:00:00 

0.1 

9 

04/21/1990  00:00:00  -  06/04/1990  23:00:00 

0.05 
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10 

07/22/1990  21:00:00  -  07/23/1990  06:00:00 

0.05 

11 

07/31/1990  13:00:00  -  07/31/1990  18:00:00 

0.05 

12 

02/05/1991  03:00:00  -  02/06/1991  23:00:00 

0.2 

13 

02/18/1991  01:00:00  -  03/22/1991  23:00:00 

0.1 

14 

03/28/1991 13:00:00  -  05/15/1991  23:00:00 

0.1 

15 

06/22/1991  07:00:00  -  09/30/1991  00:00:00 

0.05 

16 

01/08/1992  06:00:00  -  01/16/1992  23:00:00 

0.1 

17 

02/14/1992  10:00:00  -  02/16/1992  00:00:00 

0.1 
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Figure  2.  Plot  of  observed  flow  data  used  to  calibrate  the  first  HEC-HMS  model. 
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Figure  3.  Observed  flow,  quickflow,  and  base  flow  data  at  gauge  9  in  Goodwin  Creek 

Experimental  Watershed. 


Table  3  lists  the  names  and  functions  of  the  HEC-HMS  parameters  esti¬ 
mated  through  the  calibration  process.  Also  shown  in  this  table  are  the 
bounds  applied  to  these  parameters;  guidance  in  the  setting  of  these 
bounds  was  obtained  from  HEC  (2000).  In  order  to  ensure  estimation  of 
physically  acceptable  values  for  the  initial  deficit  (ID)  and  maximum  stor¬ 
age  (MS)  parameters  in  the  deficit  constant  loss  model,  the  two  adjustable 
model  parameters  LENGTH  and  A  were  related  to  the  initial  deficit  and 
maximum  storage  in  the  following  manner 

ID  =  LENGTH -A  (9) 

MS  =  LENGTH  ■  porosity  (10) 

As  noted  in  Table  3,  the  adjustable  model  parameter  A  was  specified  to  be 
less  than  the  porosity  of  the  soil,  which  was  fixed  to  be  0.36.  In  order  to 
better  accommodate  scaling  issues  resulting  from  the  use  of  different  units 
for  different  parameters,  and  in  an  attempt  to  decrease  the  degree  of 
nonlinearity  of  the  parameter  estimation  problem,  the  logs  of  these  pa¬ 
rameters  were  estimated  instead  of  their  native  values;  past  experience  has 
demonstrated  that  greater  efficiency  and  stability  of  the  parameter  estima- 
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tion  process  can  often  be  achieved  through  this  means  (Doherty  and  Ska- 
hill  2006). 


Table  3.  HEC-HMS  parameters,  their  functions,  and  constraints  imposed  during  the 

calibration  process. 


Parameter  Name 

Parameter  Function 

Bounds  Imposed  During 
Calibration 

LENGTH 

Depth  of  the  active  soil  layer 

0.003-54.681  in. 

A 

LENGTH  ■  A  =  initial  deficit;  A  <  porosity 

0.001-0.359998 

CONSTLOSS 

Infiltration  rate  when  the  soil  layer  is  satu¬ 
rated 

0.001-11.810  in./hr 

SNYDERTP 

Snyder  unit  hydrograph  standard  lag 

0.100-500  hr 

SNYDERCP 

Snyder  unit  hydrograph  peaking  coeffi¬ 
cient 

0.100  - 1 

GWSTOCOEFF 

Groundwater  storage  coefficient  for  the 
linear  reservoir 

0.100  - 10000  hr 

The  second  HEC-HMS  model  was  applied  to  the  drainage  areas  above 
streamflow  gauging  station  numbers  8,  9, 11,  and  12  in  Goodwin  Creek. 

The  second  HEC-HMS  model  included  separate  submodels  for  the  drain¬ 
age  areas  upstream  of  the  four  noted  streamflow  gauging  stations  (station 
numbers  8,  9, 11,  and  12)  located  within  the  watershed.  Precipitation  data 
from  seven  nearby  gauges  and  evaporation  data  associated  with  four  sur¬ 
rounding  locations  were  supplied  as  meteorological  forcing  data  to  sup¬ 
port  HEC-HMS  continuous  hydrologic  simulation  using  the  deficit  con¬ 
stant  loss  method. 

The  names  and  roles  of  model  parameters  selected  for  adjustment  through 
the  calibration  process  are  provided  in  Table  3.  Also  listed  are  the  bounds 
on  these  parameters  imposed  during  the  parameter  estimation  process, 
these  being  set  in  accordance  with  available  guidance  HEC  (2000).  Four 
instances  of  all  the  parameters  listed  in  Table  3  required  estimation,  one 
instance  for  each  subwatershed  model.  Thus  a  total  of  24  model  parame¬ 
ters  required  estimation  through  the  calibration  process.  As  with  the  first 
HEC-HMS  model,  in  order  to  ensure  estimation  of  physically  acceptable 
values  for  the  initial  deficit  and  maximum  storage  parameters  in  the  deficit 
constant  loss  model,  the  adjustable  model  parameters  LENGTH  and  A,  for 
each  subwatershed,  were  related  to  the  initial  deficit  and  maximum  stor¬ 
age  as  specified  in  Equations  9  and  10,  respectively.  And  as  noted  in 
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Table  3,  the  adjustable  model  parameter  A,  for  each  subwatershed,  was 
specified  to  be  less  than  the  porosity  of  the  soil,  which  was  fixed  to  be  0.36. 
Also  as  with  the  first  HEC-HMS  model,  in  order  to  better  accommodate 
scaling  issues  resulting  from  the  use  of  different  units  for  different  pa¬ 
rameters,  and  in  an  attempt  to  decrease  the  degree  of  nonlinearity  of  the 
parameter  estimation  problem,  the  logs  of  these  parameters  were  esti¬ 
mated  instead  of  their  native  values  as  past  experience  has  demonstrated 
that  greater  efficiency  and  stability  of  the  parameter  estimation  process 
can  often  be  achieved  through  this  means  (Doherty  and  Skahill  2006). 

Simultaneous  estimation  of  the  parameters  listed  in  Table  3  for  the  four 
different  subwatersheds  allows  an  important  piece  of  information  to  be 
included  in  the  parameter  estimation  process.  Namely  that,  due  to  their 
geographical  proximity  and  similarity  of  land  use,  soil  type,  and  other 
geomorphic  and  anthropogenic  conditions,  parameter  values  employed  in 
the  different  subwatershed  models  are  not  expected  to  be  significantly  dif¬ 
ferent.  To  accommodate  this  condition,  a  series  of  regularization  con¬ 
straints  effecting  an  assumed  similarity  condition  across  the  subwater¬ 
sheds  was  included  in  the  regularized  parameter  estimation  process.  That 
is,  respective  log  differences  of  identical  parameter  types  between  sub¬ 
watersheds  were  ascribed  an  “observed  value”  of  zero.  The  advantage  of 
supplying  such  information  through  regularization  constraints  rather  than 
through  “hardwired”  parameter  equality  is  that  the  regularized  inversion 
process  has  the  option  of  introducing  parameter  differences  if  this  is  a  re¬ 
quirement  for  obtaining  a  good  fit  between  modeled  and  observed  flows  at 
each  of  the  streamflow  gauging  stations.  However,  the  constrained  optimi¬ 
zation  algorithm  which  underpins  the  regularized  inversion  process  guar¬ 
antees  that  only  the  minimum  amount  of  inter-parameter  variability  re¬ 
quired  to  achieve  this  level  of  fit  is  introduced.  Thus,  subwatershed 
individuality  is  recognized  at  the  same  time  as  subwatershed  similarity. 

Estimation  of  the  24  adjustable  parameters  for  the  second  HEC-HMS 
model  was  undertaken  by  matching  observed  and  simulated  flows  over  23, 

11. 19,  and  18  non-contiguous  time  intervals  spanning  the  period  1  Jan 
1988  to  31  Dec  1990  for  stations  8,  9, 11,  and  12,  respectively.  The  noted 
periods  were  identified  based  on  a  manual  inspection  of  the  observed  flow 
data.  The  calibration  period  1  Jan  1988  to  31  Dec  1990  was  selected  based 
on  the  determination  that  there  was  no  land/use  land  cover  alteration 
within  the  four  subwatersheds  for  the  period  1988  through  1990.  Over  the 

23. 11. 19,  and  18  non-contiguous  time  intervals  spanning  the  calibration 
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period  l  Jan  1988  to  31  Dec  1990  for  stations  8,  9, 11,  and  12,  respectively, 
model  calibration  involved  matching  simulated  flows  with  observed  flows 
above  a  predetermined  threshold  (which  was  uniquely  determined  for 
each  of  the  23, 11, 19,  and  18  non-contiguous  time  intervals  for  stations  8, 
9, 11,  and  12,  respectively,  based  on  a  manual  inspection  of  the  data).  This 
resulted  in  a  total  of  8,480  flow  observations  for  use  in  the  calibration 
process  of  the  second  HEC-HMS  model. 

As  with  the  first  HEC-HMS  model,  the  objective  function  for  the  second 
HEC-HMS  model  was  defined  as  the  sum  of  weighted  squared  differences 
between  modeled  and  observed  log-transformed  flows,  with  all  weights 
assigned  a  value  of  1.0.  Thus  h  of  Equation  3  was  comprised  of  the  logs  of 
flows,  while  the  model  represented  by  X  in  these  equations  calculated  the 
model -generated  counterparts  to  these  logged  flows.  Q  was  the  identity 
matrix.  For  the  second  HEC-HMS  model,  the  simulation  time  interval  was 
one  hour,  which  equaled  the  temporal  resolution  of  the  input  precipitation 
data. 

Tables  4  to  7  list  the  23, 11, 19,  and  18  non-contiguous  time  intervals  for 
stations  8,  9, 11,  and  12,  respectively,  and  their  respective  thresholds.  Fig¬ 
ures  4  to  7  are  plots  of  the  observed  flow  data  that  were  compared  with 
simulated  flows  to  calibrate  the  second  HEC-HMS  model. 


Table  4.  For  the  second  HEC-HMS  model,  non-contiguous  time  intervals  and  thresholds 

associated  with  station  8. 


Number 

Time  Interval 

Threshold  (cfs) 

1 

01/18/1988  20:00:00  -  01/20/1988  00:00:00 

0.1 

2 

02/02/1988  02:00:00  -  02/05/1988  00:00:00 

0.1 

3 

02/14/1988  16:00:00  -  03/14/1988  00:00:00 

0.2 

4 

04/06/1988  00:00:00  -  04/07/1988  00:00:00 

0.1 

5 

11/19/1988  20:00:00  - 11/27/1988  00:00:00 

0.2 

6 

12/21/1988  00:00:00  -  01/02/1989  00:00:00 

0.2 

7 

01/07/1989  16:00:00  -  02/04/1989  10:00:00 

0.1 

8 

02/12/1989  18:00:00  -  03/08/1989  00:00:00 

0.1 

9 

03/29/1989  04:00:00  -  04/06/1989  00:00:00 

0.1 

10 

05/05/1989  02:00:00  -  05/23/1989  00:00:00 

0.3 

11 

06/08/1989  02:00:00  -  06/17/1989  00:00:00 

0.2 

12 

07/01/1989  06:00:00  -  07/17/1989  00:00:00 

0.2 

13 

07/30/1989  20:00:00  -  08/07/1989  16:00:00 

0.1 

14 

11/06/1989  00:00:00  - 11/09/1989  00:00:00 

0.2 
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15 

12/30/1989  02:00:00  -  01/08/1990  12:00:00 

0.2 

16 

01/17/1990  14:00:00  -  01/30/1990  00:00:00 

0.3 

17 

02/02/1990  00:00:00  -  03/18/1990  00:00:00 

0.2 

18 

03/28/1990  12:00:00  -  04/08/1990  00:00:00 

0.2 

19 

04/21/1990  08:00:00  -  05/22/1990  12:00:00 

0.3 

20 

06/03/1990  00:00:00  -  06/04/1990  00:00:00 

0.2 

21 

11/28/1990  00:00:00  - 12/04/1990  00:00:00 

0.2 

22 

12/17/1990  00:00:00  - 12/23/1990  00:00:00 

0.2 

23 

12/27/1990  04:00:00  - 12/31/1990  23:00:00 

0.2 

Table  5.  For  the  second  HEC-HMS  model,  non-contiguous  time  intervals  and  thresholds 

associated  with  station  9. 


Number 

Time  Interval 

Threshold  (cfs) 

1 

01/02/1989  00:00:00  -  02/13/1989  23:00:00 

0.1 

2 

02/18/1989  06:00:00  -  02/21/1989  12:00:00 

0.1 

3 

02/27/1989  02:00:00  -  05/17/1989  23:00:00 

0.1 

4 

05/22/1989  00:00:00  -  08/08/1989  23:00:00 

0.1 

5 

01/03/1990  18:00:00  -  01/08/1990  08:00:00 

0.1 

6 

01/28/1990  20:00:00  -  01/29/1990  23:00:00 

0.1 

7 

02/09/1990  13:00:00  -  03/13/1990  23:00:00 

0.1 

8 

04/06/1990  00:00:00  -  04/06/1990  23:00:00 

0.1 

9 

04/21/1990  00:00:00  -  06/04/1990  23:00:00 

0.05 

10 

07/22/1990  21:00:00  -  07/23/1990  06:00:00 

0.05 

11 

07/31/1990  13:00:00  -  07/31/1990  18:00:00 

0.05 

Table  6.  For  the  second  HEC-HMS  model,  non-contiguous  time  intervals  and  thresholds 

associated  with  station  11. 


Number 

Time  Interval 

Threshold  (cfs) 

1 

01/18/1988  20:00:00  -  01/21/1988  00:00:00 

0.1 

2 

02/02/1988  02:00:00  -  03/14/1988  00:00:00 

0.1 

3 

04/02/1988  00:00:00  -  04/13/1988  00:00:00 

0.1 

4 

11/26/1988  02:00:00  -  02/04/1989  00:00:00 

0.3 

5 

02/12/1989  18:00:00  -  03/06/1989  10:00:00 

0.2 

6 

03/29/1989  04:00:00  -  04/05/1989  12:00:00 

0.1 

7 

05/05/1989  02:00:00  -  05/11/1989  00:00:00 

0.1 

8 

06/05/1989  16:00:00  -  06/17/1989  00:00:00 

0.1 

9 

07/01/1989  04:00:00  -  07/17/1989  00:00:00 

0.1 
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10 

07/30/1989  20:00:00  -  08/08/1989  00:00:00 

0.1 

11 

11/06/1989  00:00:00  - 11/09/1989  00:00:00 

0.1 

12 

12/30/1989  00:00:00  -  01/08/1990  10:00:00 

0.2 

13 

01/28/1990  18:00:00  -  03/03/1990  10:00:00 

0.2 

14 

03/07/1990  12:00:00  -  03/17/1990  00:00:00 

0.3 

15 

03/28/1990  10:00:00  -  04/07/1990  12:00:00 

0.2 

16 

04/21/1990  06:00:00  -  05/22/1990  10:00:00 

0.2 

17 

06/03/1990  02:00:00  -  06/04/1990  12:00:00 

0.1 

18 

11/27/1990  22:00:00  - 12/04/1990  00:00:00 

0.1 

19 

12/17/1990  00:00:00  - 12/31/1990  23:00:00 

0.3 

Table  7.  For  the  second  HEC-HMS  model,  non-contiguous  time  intervals  and  thresholds 

associated  with  station  12. 


Number 

Time  Interval 

Threshold  (cfs) 

1 

01/17/1988  00:00:00  -  01/20/1988  12:00:00 

0.2 

2 

02/02/1988  00:00:00  -  02/05/1988  00:00:00 

0.1 

3 

02/14/1988  18:00:00  -  04/12/1988  12:00:00 

0.1 

4 

11/19/1988  20:00:00  -  11/27/1988  12:00:00 

0.1 

5 

12/21/1988  00:00:00  -  12/31/1988  23:00:00 

0.1 

6 

01/11/1989  01:00:00  -  02/28/1989  23:00:00 

0.1 

7 

05/05/1989  01:00:00  -  05/10/1989  12:00:00 

0.2 

8 

06/04/1989  06:00:00  -  06/17/1989  00:00:00 

0.1 

9 

07/01/1989  04:00:00  -  07/17/1989  12:00:00 

0.1 

10 

08/06/1989  16:00:00  -  08/08/1989  00:00:00 

0.1 

11 

11/06/1989  00:00:00  -  11/09/1989  00:00:00 

0.1 

12 

12/30/1989  00:00:00  -  01/08/1990  12:00:00 

0.2 

13 

01/28/1990  18:00:00  -  02/04/1990  00:00:00 

0.4 

14 

02/09/1990  12:00:00  -  03/17/1990  00:00:00 

0.3 

15 

03/28/1990  08:00:00  -  04/07/1990  00:00:00 

0.3 

16 

04/21/1990  08:00:00  -  06/10/1990  00:00:00 

0.2 

17 

11/27/1990  22:00:00  -  11/28/1990  12:00:00 

0.1 

18 

12/17/1990  02:00:00  - 12/31/1990  23:00:00 

0.2 
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Figure  4.  Plot  of  station  8  observed  flow  data  used  to  calibrate  the  second  HMS  model. 
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Figure  5.  Plot  of  station  9  observed  flow  data  used  to  calibrate  the  second  HMS  model 
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Figure  6.  Plot  of  station  11  observed  flow  data  used  to  calibrate  the  second  HMS  model. 


Figure  7.  Plot  of  station  12  observed  flow  data  used  to  calibrate  the  second  HMS  model. 
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5  Results 

First  HEC-HMS  Model  -  Calibration  Experiment  1 

In  an  attempt  to  locate  as  many  local  minima  as  possible,  PD_MS2  was 
asked  to  run  loo  inversion  runs  from  a  succession  of  starting  values  that 
were  maximally  distant  from  all  previous  parameter  trajectories,  as  dis¬ 
cussed  above.  These  starting  values  were  selected  from  1,296  random  pa¬ 
rameter  samples  for  which  objective  functions  were  calculated  prior  to  the 
undertaking  of  any  inversion  runs. 

Figure  8  depicts  the  outcomes  of  this  exercise.  Each  of  the  six  graphs  ap¬ 
pearing  in  this  figure  pertains  to  one  of  the  six  estimated  parameters.  For 
each  graph  the  objective  function  is  plotted  on  the  x  axis,  while  an  opti¬ 
mized  parameter  value  is  plotted  on  the  y  axis.  In  each  graph,  each  point 
represents  the  outcome  of  one  PD_MS2-supervised  inversion  run;  corre¬ 
sponding  points  from  different  graphs  (representing  corresponding  values 
for  different  parameters)  can  be  matched  vertically  through  their  common 
objective  function.  It  is  readily  apparent  from  this  figure  that  many  of  the 
outcomes  of  successive  optimization  runs  are  grouped  into  “parameter 
clumps”  of  nearly  constant  objective  function  value;  these  clumps  define 
regions  of  attraction  in  parameter  space.  “Tight”  clumps  indicate  a  well- 
defined  region  of  attraction;  vertical  spreading  of  clumps  indicates  diffi¬ 
culties  in  parameter  identification  through  parameter  correlation  and/or 
insensitivity.  For  those  minima  situated  at  the  bottom  of  broad  objective 
function  valleys  defining  different  regions  of  attraction  in  parameter  space, 
local  minima  are  often  in  close  proximity.  Other  local  minima  appear  to 
exist  in  isolation  from  these  more  populous  clumps. 

For  this  calibration  problem,  the  objective  function  has  a  value  of  223.6  at 
its  global  minimum.  Table  8  summarizes  computed  Nash-Sutcliffe  effi¬ 
ciency  scores  (which  are  based  on  a  comparison  of  the  flow  observations 
used  in  the  calibration  process)  for  each  of  the  individual  17  non¬ 
contiguous  time  intervals.  Figure  9  shows  the  fit  between  modeled  and  ob¬ 
served  flows  at  this  minimum  for  the  flows  that  were  compared  for  the 
seventh  non-contiguous  time  interval  listed  in  Table  2. 
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Figure  8.  End-points  in  parameter  space  of  100  inversion  runs  undertaken  under  the  control 
of  PD_MS2.  Parameters  comprising  an  optimized  set  are  linked  vertically  between  graphs  by 
objective  function.  TPI  functionality  was  not  operative. 


Ninety-five  percent  confidence  intervals  for  the  parameters  corresponding 
to  the  global  objective  function  minimum  are  provided  in  Table  9,  and 
these  could  only  be  calculated  after  fixing  the  parameter  A,  for  A  could  not 
be  defined  at  the  global  objective  function  minimum.  The  insensitivity  of 
parameter  A  is  exhibited  upon  inspection  of  Figure  8,  wherein  it  is  appar¬ 
ent  that  the  global  objective  function  minimum  of  223.6  occurs  in  an  elon¬ 
gate  valley  rather  than  a  bowl  (see  the  dotty  plot  for  the  initial  deficit  pa¬ 
rameter,  ID). 
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Figure  9.  Simulated  and  observed  flows  at  identified  global  minimum  for  the  flows  that  were 
compared  for  the  seventh  non-contiguoustime  interval  listed  in  Table  2. 


Table  8.  Nash-Sutcliffe  efficiency  associated  with  identified  global  minimum. 


Number 

Number  of  Series  Terms  in  this  Interval 

Nash-Sutcliffe  Coefficient 

1 

175 

0.61 

2 

34 

0.32 

3 

221 

0.18 

4 

196 

0.24 

5 

81 

0.60 

6 

25 

0.83 

7 

210 

0.72 

8 

24 

0.89 

9 

155 

0.38 

10 

5 

0.61 

11 

4 

0.13 

12 

32 

0.71 

13 

248 

0.51 

14 

463 

0.44 

15 

61 

0.22 

16 

90 

0.59 

17 

29 

-0.23 
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Table  9.  Parameter  values  corresponding  to  global  optimum;  also  shown  are  linear 
parameter  confidence  limits  calculated  as  a  by-product  of  the  GML  parameter  estimation 

process. 


Parameter  name 

Estimated  Value 

Lower  95  Percent 
Confidence  Limit 

Upper  95  Percent 
Confidence  Limit 

LENGTH 

0.816390 

0.793602 

0.839832 

CONSTLOSS 

0.151830 

0.143696 

0.160425 

SNYDERTP 

126.100 

123.849 

128.392 

SNYDERCP 

0.100000 

9.711740E-02 

0.102968 

GWSTOCOEFF 

1.87410 

1.79730 

1.95418 

Further  evidence  supporting  the  existence  of  local  optima  for  the  first 
HEC-HMS  model  was  provided  through  optimization  runs  performed 
within  HEC-HMS  using  the  Nelder  and  Mead  (1965)  local  search  method 
and  the  sum  of  squared  residuals  objective  function.  Table  10  summarizes 
the  results  of  the  optimization  runs. 


Table  10.  Initial  parameter  sets,  final  parameter  sets,  and  objective  function  values 
associated  with  three  HEC-HMS  optimization  runs  with  the  first  HEC-HMS  model. 


!  TRIAL  I 

1  1 

1  2  1 

1  3 

Parameter  List 

Initial  Deficit 

Max  Storage 
Constant  Loss 
Snyder  Peak 
Snyder  Coeff 
GW  Coeff 

Initial  Parameter 

Value 

0.0063 

0.0137 

6.9169 

385.27 

0.4062 

1165.4 

HMS  Optimized 
Parameter  Value 

0.1032 

0.4642 

7.0139 

385.37 

0.5032 

1165.5 

PD_MS2 
Optimized 
Parameter  Value 

0.0005 

0.0010 

1 .6772 

500 

0.6882 

853.75 

Initial  Parameter 

Value 

3.9775 

5.0468 

2.7432 

0  2061 

0.7416 

9470.9 

HMS  Optimized 
Parameter  Value 

3.9775 

5.0468 

2.7432 

0.2061 

0.7416 

9470.9 

PD_MS2 
Optimized 
Parameter  Value 

0.0134 

0.0918 

0.0416 

0  2061 

0.7416 

11.83 

Initial  Parameter 

Value 

0.0667 

0  5799 

8.5194 

298.48 

0.1470 

9738.8 

HMS  Optimized 
Parameter  Value 

0.1636 

1.0304 

8.6164 

298.58 

0  2440 
9738.9 

PD_MS2 
Optimized 
Parameter  Value 

0.0130 

0.5673 

8.5194 

298.48 

0.1470 

8083.3 

HMS  Objective 
Function  Value 

49691.1000 

50275.7000 

50170.1000 

Figure  10  is  a  plot  of  simulated  flows  and  base  flows  at  the  global  objective 
function  minimum.  While  we  are  now  able  to  identify  the  global  objective 
function  minimum  using  automated  parameter  estimation  capabilities 
with  a  HEC-HMS  model,  this  new  capability  does  not  preclude  identifica¬ 
tion  of  a  hydrologically  unacceptable  model  at  the  global  minimum.  A 
more  robust  objective  function  formulation  is  needed  that  better  reflects 
our  perceptual  model  of  the  watershed  system. 
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First  HEC-HMS  Model  -  Calibration  Experiment  2 

In  an  attempt  to  locate  as  many  local  minima  as  possible,  PD_MS2  was  asked 
to  run  100  inversion  runs  from  a  succession  of  starting  values  that  were  maxi¬ 
mally  distant  from  all  previous  parameter  trajectories,  as  discussed  above. 
These  starting  values  were  selected  from  1,296  random  parameter  samples  for 
which  objective  functions  were  calculated  prior  to  the  undertaking  of  any  in¬ 
version  runs.  Figure  11  depicts  the  outcomes  of  this  exercise. 

For  this  calibration  problem,  the  objective  function  has  a  value  of  428.6  at 
its  global  minimum.  Table  11  summarizes  computed  Nash-Sutcliffe  effi¬ 
ciency  scores  (which  are  based  on  a  comparison  of  the  flow  observations 
used  in  the  calibration  process)  for  each  of  the  individual  17  non¬ 
contiguous  time  intervals.  Figure  12  shows  the  fit  between  modeled  and 
observed  flows  at  this  minimum  for  the  flows  that  were  compared  for  the 
seventh  non-contiguous  time  interval  listed  in  Table  2. 

Figure  13  is  a  plot  of  the  simulated  flows  and  base  flows  at  the  global  ob¬ 
jective  function  minimum.  While  we  have  now  demonstrated  how  state 
information  other  than  stream  discharge  data  alone  can  be  included  into 
the  automatic  calibration  process  of  a  HEC-HMS  model,  which  results,  at 
least  for  this  case,  in  a  much  more  hydrologically  acceptable  model  than 
calibrating  solely  against  flows,  this  process  still  may  not  remove  the  prob¬ 
lem  of  local  optima. 
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Figure  11.  End-points  in  parameter  space  of  100  inversion  runs  undertaken  under  the 
control  of  PD_MS2.  Parameters  comprising  an  optimized  set  are  linked  vertically  between 
graphs  by  objective  function.  TPI  functionality  was  not  operative. 
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Table  11.  Nash-Sutcliffe  efficiency  associated  with  identified  global  minimum. 


Number 

Number  of  Series  Terms  in  this  Interval 

Nash-Sutcliffe  Coefficient 

1 

175 

0.61 

2 

34 

0.70 

3 

221 

0.55 

4 

196 

0.51 

5 

81 

0.76 

6 

25 

0.09 

7 

210 

0.70 

8 

24 

0.73 

9 

155 

0.56 

10 

5 

-14.18 

11 

4 

-9.73 

12 

32 

0.34 

13 

248 

0.73 

14 

463 

0.74 

15 

61 

-0.13 

16 

90 

0.33 

17 

29 

-0.20 
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Figure  12.  Simulated  and  observed  flows  at  identified  global  minimum  for  the  flows  that  were 
compared  for  the  seventh  non-contiguoustime  interval  listed  in  Table  2. 
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Second  HEC-HMS  Model 

Table  12  lists  parameter  values  for  each  subwatershed  model,  estimated 
using  the  adaptive  regularization  scheme  described  above.  In  implement¬ 
ing  the  regularized  inversion  process,  a  very  low  target  measurement  ob¬ 
jective  function  was  set;  hence  Om  of  Equation  3  was  lowered  as  far  as  pos¬ 
sible,  thus  reducing  misfit  between  measured  and  observed  flows  to  a 
minimum.  It  is  apparent  from  Table  12  that  optimal  fitting  of  model  out¬ 
puts  to  matched  flows  could  only  be  achieved  through  the  assignment  of 
different  values  to  parameters  of  the  same  type  in  different  sub  watersheds. 
However,  the  adaptive  regularization  scheme  employed  in  their  estimation 
attempted  to  ensure  that  these  differences  were  kept  to  a  minimum.  The 
total  measurement  objective  function  (pertaining  to  all  streamflow  gauges) 
achieved  through  this  calibration  exercise  was  1,047. 

Table  13  summarizes  computed  Nash-Sutcliffe  efficiency  scores  (which  are 
based  on  a  comparison  of  the  flow  observations  used  in  the  calibration 
process)  for  each  of  the  individual  non-contiguous  time  intervals  for  each 
gauged  subwatershed. 

Figure  14  shows  the  fit  between  the  modeled  and  observed  flows  for 
streamflow  gauging  station  9  for  the  flows  that  were  compared  for  the 
seventh  non-contiguous  time  interval  listed  in  Table  2. 
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Table  12.  Estimated  values  for  subwatershed  model  parameters  for  attainment  of  best  fit  at  all 
subwatershed  streamflow  gauging  stations,  this  corresponding  to  a  measurement  objective 
function  of  135.1.  Adaptive  regularization  was  employed  in  the  parameter  estimation  process. 


Parameter  Name 

Station  8 

Station  9 

Station  11 

Station  12 

MS 

1.031548 

0.193263 

0.169596 

0.209377 

ID 

0.143271 

0.026842 

0.023555 

0.029080 

CONSTLOSS 

0.024030 

0.039661 

0.080199 

0.039745 

SNYDERTP 

0.250000 

0.250000 

0.250000 

0.250000 

SNYDERCP 

0.560000 

0.560000 

0.560000 

0.560000 

GWSTOCOEFF 

12.316100 

16.01440 

11.648300 

9.545690 

Table  13.  Number  of  terms  and  Nash-Sutcliffe  efficiency  scores,  ES,  associated  with  identified 

minimum  for  each  station. 


Number 

Station 

8 

9 

11 

12 

Measure¬ 

ment 

ES 

Measure¬ 

ment 

ES 

Measure¬ 

ment 

ES 

Measure¬ 

ment 

ES 

1 

29 

0.76 

175 

0.79 

45 

0.85 

58 

0.74 

2 

69 

-2.35 

34 

0.82 

306 

0.17 

306 

0.17 

3 

237 

0.53 

221 

0.69 

84 

0.75 

84 

0.75 

4 

18 

0.35 

196 

0.69 

215 

-0.04 

356 

0.12 

5 

25 

-18.92 

81 

0.69 

249 

0.86 

334 

0.86 

6 

95 

-0.28 

25 

0.50 

104 

0.64 

104 

0.64 

7 

288 

0.66 

210 

0.82 

63 

0.68 

58 

0.68 

8 

393 

0.81 

24 

0.66 

123 

0.46 

123 

0.46 

9 

124 

0.28 

155 

0.62 

94 

0.36 

94 

0.36 

10 

115 

0.81 

5 

-20.57 

45 

-2.66 

45 

-2.66 

11 

101 

0.66 

4 

-12.35 

27 

-13.11 

27 

-13.11 

12 

26 

0.75 

132 

0.63 

132 

0.63 

13 

69 

-2.35 

208 

0.86 

142 

0.86 

14 

237 

0.53 

95 

0.72 

95 

0.72 

15 

14 

0.27 

92 

-1.52 

65 

-2.07 

16 

24 

-19.59 

111 

-0.30 

111 

-0.30 

17 

95 

-0.28 

29 

-12.04 

29 

-12.04 

18 

248 

0.66 

20 

-4.26 

14 

-5.09 

19 

327 

0.81 

158 

0.52 

20 

108 

0.26 

45 

21 

129 

0.81 

306 

22 

101 

0.66 

84 

23 

101 

0.66 

215 
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Figure  14.  Simulated  and  observed  flows  at  station  9,  at  the  minimum  objective  function 
value,  for  the  flows  that  were  compared  for  the  seventh  non-contiguous  time  interval  listed 

in  Table  2. 

Figure  15  is  a  plot  of  the  simulated  flows  and  base  flows  for  streamflow 
gauging  station  9  at  the  objective  function  minimum. 


Figure  15.  Plot  of  simulated  flows,  quickflow,  and  base  flow,  for  station  9,  at  the  minimum. 
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6  Discussion 

The  objectives  for  this  article  were  to  describe  and  demonstrate  the  use  of 
parameter  estimation  methodologies  that  could  potentially  be  employed  to 
improve  upon  existing  HEC-HMS  (Hydrologic  Engineering  Center’s  Hydro- 
logic  Modeling  System)  automated  parameter  estimation  capabilities.  In 
particular,  the  intent  for  this  article  was  to  describe  and  demonstrate  the 
use  of  methods  that  (l)  accommodate  local  minima  (Skahill  and  Doherty 
2006)  and  (2)  support  the  inversion  of  complex  (i.e.,  highly  parameter¬ 
ized)  models  (Doherty  and  Skahill  2006),  and  through  this  process  to  also: 

a.  Demonstrate  that  each  of  these  two  methods  provide  information 
about  individual  parameter  sensitivities  and  parameter  correlation. 

b.  Demonstrate  how  state  information  other  than  stream  discharge  data 
alone  can  be  included  into  the  automatic  calibration  process  of  an 
HEC-HMS  model. 

c.  Demonstrate  how  objective  function  definitions  associated  with  spe¬ 
cific  periods  of  the  calibration  dataset  can  be  included  into  the  auto¬ 
matic  calibration  process  of  an  HEC-HMS  model. 

The  principal  intent  of  the  first  calibration  experiment  with  the  first  HEC- 
HMS  model  was  to  illustrate  the  existence  of  multiple  local  optima  (and 
this  was  also  demonstrated  through  the  trials  using  the  automatic  calibra¬ 
tion  capabilities  currently  available  within  HEC-HMS),  and  to  also  dem¬ 
onstrate  that  the  trajectory  repulsion  scheme  can  be  used  with  HEC-HMS 
to  accommodate  the  presence  of  multiple  local  optima  and  find  the  global 
objective  function  minimum. 

The  initial  deficit  parameter  of  the  deficit  constant  loss  model  (HEC  2005) 
was  identified  to  be  completely  insensitive  at  the  identified  global  objective 
function  minimum  for  the  first  HEC-HMS  model.  Estimates  of  parameter 
uncertainty,  correlation  and  (in)sensitivity  are  readily  available  as  a  by¬ 
product  of  the  use,  both  during  and  at  the  end  of  a  model  inversion,  of  the 
parameter  estimation  methods  described  in  the  theory  section  of  this 
article,  and  such  information  can  be  readily  employed  to  reformulate  a 
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problem  to  achieve  a  more  stable  inversion  and  likely  also  a  more  optimal 
model  estimate. 

Simulated  flows  from  the  estimated  model  at  the  global  objective  function 
minimum  for  the  first  calibration  experiment  with  the  first  HEC-HMS 
model  fit  observed  stream  discharges  fairly  well;  however,  upon  further 
inspection  it  was  an  unacceptable  model  in  that  almost  all  simulated  flow 
was  base  flow  (see  Figure  10).  This  contradicted  our  perceptual  model  for 
the  watershed  system,  and  past  research,  wherein  system  response  is 
dominated  by  direct  runoff.  This  observation  underscores  the  fact  that  the 
ability  to  find  the  global  objective  function  minimum  is  an  insufficient  re¬ 
quirement  to  attain  a  hydrologically  acceptable  model.  If  at  all  possible, 
one  must  process  the  observation  dataset(s)  in  a  manner  that  best  reflects 
our  perceptual  understanding  of  a  system’s  hydrologic  response. 

The  principal  intent  of  the  second  calibration  experiment  with  the  first 
HEC-HMS  model  was  to  demonstrate  that  state  information  other  than 
stream  discharge  data  alone  can  be  included  into  the  automatic  calibration 
process  of  an  HEC-HMS  model,  and  through  that  process,  at  least  for  this 
case,  one  can  attain  a  more  hydrologically  acceptable  model  than  calibrat¬ 
ing  solely  against  flows  (see  Figure  13),  and  moreover,  retain  a  reasonable 
fit  with  the  actual  hard  data  (see  Tables  8  and  11  and  Figures  9  and  12). 
However,  as  evidenced  by  inspecting  Figure  11,  this  process  does  not  cir¬ 
cumvent  the  presence  of  multiple  local  optima. 

The  calibration  experiment  with  the  second  HEC-HMS  model  demon¬ 
strated  that  one  can  simultaneously  calibrate  multiple  subwatershed  mod¬ 
els  automatically  with  an  HEC-HMS  model.  The  strength  of  such  an  ap¬ 
proach  is  that  it  includes  into  the  calibration  process  explicit  recognition  of 
the  fact  that  the  semi-physical  basis  of  parameters  employed  by  models 
such  as  HEC-HMS  demands  that  parameter  values  associated  with  similar 
land  uses  and  soil  types  in  adjacent  areas  be  at  least  broadly  similar.  Man¬ 
ual  calibration  of  the  four  Goodwin  Creek  subwatershed  models,  if  care¬ 
fully  implemented,  would  probably  follow  this  strategy  (i.e.,  calibrate  each 
model  individually,  with  due  recognition  of  the  desirability  of  inter¬ 
subwatershed  parameter  similarity),  for  a  modeler  would  be  aware  of  the 
need  for  maintaining  at  least  a  certain  degree  of  parameter  similarity 
across  subwatersheds  as  he/she  assigns  values  to  them  in  successive  trial- 
and-error  model  runs.  However,  it  is  likely  that  this  would  be  a  difficult 
undertaking  from  a  practical  point  of  view,  probably  requiring  much  pa- 
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tience  on  the  part  of  the  modeler,  while  many  model  runs  are  undertaken 
to  achieve  this  ideal.  Furthermore  such  an  exercise  would  leave  the  mod¬ 
eler  with  no  knowledge  of  whether  model -to-measurement  fits  on  the  one 
hand,  and  inter-subwatershed  parameter  similarity  on  the  other  hand, 
could  be  further  improved  with  an  even  greater  amount  of  time  devoted  to 
the  manual  calibration  exercise.  Hence,  computer  assistance  in  this  proc¬ 
ess  is  obviously  desirable.  However,  the  use  of  automated  methods  of 
model  calibration  must  not  reduce  the  capacity  of  the  modeler  to  exercise 
his/her  judgment  in  this  and  other  matters;  in  fact,  it  should  enhance  it. 
Fortunately,  regularized  inversion  does  indeed  allow  a  modeler’s  judgment 
to  become  an  integral  part  of  the  calibration  process  through  implement¬ 
ing  that  process  as  a  constrained  minimization  problem,  with  the  con¬ 
straints  being  the  modeler’s  idea  of  what  constitutes  an  ideal  parameter 
set.  Departures  from  these  constraints  are  tolerable  only  to  the  extent  that 
they  result  in  an  “adequate”  (as  defined  by  the  modeler)  level  of  model-to- 
measurement  fit. 

For  the  Goodwin  Creek  Experimental  Watershed,  inclusion  of  a  modeler’s 
wisdom  in  the  calibration  process  possibly  assumes  more  importance  than 
in  many  other  modeling  contexts  because  of  the  potential  paucity  of  data 
available  for  the  calibration  of  each  individual  subwatershed,  and  hence 
the  high  potential  for  parameter  nonuniqueness.  Simultaneous  calibration 
of  four  subwatershed  models  allows  the  quantity  of  data  employed  by  the 
parameter  estimation  process  to  be  increased  by  a  factor  of  four.  In  theory, 
this  provides  some  relief  from  the  effects  of  data  paucity  in  estimating  pa¬ 
rameters  for  four  individual  subwatershed  models  through  separate  cali¬ 
bration  exercises.  This  is  because  it  does  this  without  necessarily  introduc¬ 
ing  four  times  the  number  of  parameters  to  the  calibration  process,  for  the 
constrained  optimization  process  implemented  through  regularized  inver¬ 
sion  dispenses  with  the  need  for  parameters  in  different  subwatersheds  to 
vary  independently  of  each  other  unless  the  information  content  of  the 
calibration  dataset  makes  such  parameter  differences  tenable.  Where  a 
certain  degree  of  parameter  independence  is  proven  to  be  warranted,  extra 
“observations”  are  introduced  to  the  parameter  estimation  process 
through  the  regularization  algorithm  through  which  this  process  is  imple¬ 
mented  -  “observations”  whose  presence  should  endow  that  process  with 
extra  numerical  stability  because  of  the  fact  that  they  pertain  directly  to 
model  parameters  themselves,  and  that  at  least  one  such  “observation”  ex¬ 
ists  for  every  estimable  parameter.  In  theory,  the  result  is  a  stable  process 
that  allows  maximum  receptivity  of  parameters  to  both  “hard  information” 
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provided  by  the  measurement  dataset  and  “soft  information”  embodied  in 
a  modeler’s  understanding  of  the  area,  encapsulated  in  the  set  of  regulari¬ 
zation  constraints.  The  results  shown  in  Table  13  (in  comparison  to  Tables 
8  and  11)  and  Figures  14  and  15  seem  to  confirm  the  above  noted  com¬ 
ments. 

All  three  calibration  experiments  demonstrated  how  objective  function 
definitions  associated  with  specific  periods  of  the  calibration  dataset  can 
be  included  into  the  automatic  calibration  process  of  an  HEC-HMS  model, 
and  this  was  achieved  due  to  the  model  independent  nature  of  the  meth¬ 
ods  described. 

While  the  PD_MS2  “multi-start”  procedure  is  smart  in  that  it  uses  the  tra¬ 
jectories  from  previous  local  searches  to  select  new  initial  conditions  that 
are  maximally  distant  in  parameter  space,  it  is  nevertheless  inefficient  in 
that  it  explores  many  parts  of  parameter  space  which  are  dead  ends  (and 
this  is  clearly  evident  upon  inspection  of  Figures  8  and  11).  Hence,  future 
work  will  focus  on  the  selection,  adaptation,  and  implementation  of  a  hy¬ 
brid  genetic  and  gradient-based  optimization  algorithm  for  identifying 
globally  optimal  model  parameters. 

Further  research  is  also  needed  to  address  proper  objective  function  for¬ 
mulation,  and  weights  assignment,  to  ensure  that  the  model  identified  at 
the  conclusion  of  an  estimation  process  is  reflective  of  the  multiple  proc¬ 
esses  (which  operate  at  different  time  scales)  encapsulated  in  the  model 
structure,  intended  uses  for  the  model,  and  the  data  available  to  identify 
the  parameters  related  to  those  processes  and  intended  predictions.  An¬ 
other  avenue  for  investigation  will  be  the  use  of  multi-resolution  analysis 
(wavelets)  to  decompose  the  observed  system  response  data  into  different 
time  scales  to  reflect  the  time  scales  of  differing  processes  in  the  watershed 
model. 

Further  research  will  also  explore  the  use  of  additional  regularization 
strategies.  Utilization  of  the  time-scale  decomposed  objective  function  will 
allow  regularization  strategies  that  reflect  the  physical  time  scales  present 
in  the  model  and  observed  response  data. 
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